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The  Finite  Difference  Time  Domain  (FDTD)  method  is  limited  by  memory 
requirements  and  computation  time  when  applied  to  large  problems,  complicated 
geometries,  or  geometries  with  fine  features.  In  this  thesis,  the  nonuniform  orthogonal 
FDTD  method  is  presented  and  applied  to  a  variety  of  electromagnetic  problems.  The 
nonuniform  aspect  of  the  method  gives  great  flexibility  in  modeling  complicated 
geometries  with  fine  features.  Furthermore,  the  variability  of  the  mesh  resolution  also 
enables  the  user  to  move  the  boundaries  of  the  computational  domain  farther  away  from 
the  center  of  the  problem  without  an  undue  increase  in  the  number  of  cells.  Most 
significantly,  the  orthogonality  of  the  method  preserves  the  speed  of  the  conventional 
FDTD  method.  These  three  features  of  the  nonuniform  orthogonal  FDTD  method  are 
demonstrated  by  means  of  numerical  examples  throughout  the  thesis. 

Grid  dispersion  error  from  the  nonuniform  mesh  is  analyzed  and  numerical 
examples  are  presented,  demonstrating  that  small  growth  rates  in  mesh  discretization  lead 
to  acceptably  small  errors.  The  issue  of  absorbing  boundary  conditions  is  addressed  with 
the  analysis  and  application  of  the  dispersive  boundary  condition  on  nonuniform  meshes. 
New  techniques  are  also  introduced  for  the  efficient  characterization  of  microstrip  lines, 
microstrip  discontinuities,  and  coupled  microstrip  structures  using  FDTD  data.  A  local 
mesh  refinement  technique  is  introduced  for  planar  perfect  electric  conductor,  and  is 
shown  to  be  three  times  more  accurate  than  the  staircasing  approximation. 

The  versatility  of  the  method  is  demonstrated  by  the  analysis  of  a  balun-fed  folded 
dipole  antenna,  the  characterization  of  the  transition  of  grounded  coplanar  waveguide  to 
microstrip  line,  and  the  study  of  fields  in  lossy  layered  media. 
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CHAPTER  1 


INTRODUCTION 


1.1  Overview 

The  Finite  Difference  Time  Domain  (FDTD)  method  was  introduced  to  the 
electromagnetics  community  in  1966  [1.1].  The  FDTD  method  solves  Maxwell's 
equations  in  partial  differential  equation  form.  It  is  extremely  versatile  in  that  it  is  able  to 
analyze  geometrical  structures  with  arbitrary  inhomogeneities.  This  is  an  attractive 
aspect  of  the  method  when  compared  with  integral  equation  methods,  which  require 
geometry-dependent  Green's  functions;  these  can  be  computationally  expensive  and 
difficult  to  derive  for  arbitrary  inhomogeneities.  Although  approximations  often  exist,  it 
frequently  remains  difficult  to  obtain  robust  representations  of  the  Green's  functions  for 
the  general  case  [1.2].  Moreover,  the  FDTD  method  does  not  require  matrix  operations, 
unlike  the  Finite  Element  Method  and  the  Method  of  Moments.  Depending  upon  the 
problem,  the  matrix  operations  present  potential  pitfalls  for  those  methods  because  of 
ill-conditioned  matrices,  non-convergent  solutions,  or  simply  the  computational  expense 
of  N  by  N  matrix  operations.  Additionally,  the  FDTD  method  is  inherently  efficient 
since  it  is  of  order  N  [1.3]. 

The  FDTD  method  is  explicit  and  solves  differential  equations  by  stepping 
forward  in  time.  Because  it  is  a  time  domain  method,  solutions  to  transient  excitations 
provide  a  wide  band  of  results  in  the  frequency  domain.  This  is  particularly  useful  when 
analyzing  frequency-dependent  structures,  since  the  frequency-dependent  results  can  be 
obtained  from  a  single  FDTD  simulation.  Furthermore,  the  FDTD  method  is  an  excellent 
tool  for  providing  a  means  for  field  visualization  as  a  function  of  time  and  space.  This 
can  be  important  when  studying  complex  three-dimensional  geometries,  where  analysis 
of  time  signatures  of  fields,  voltages,  or  currents  can  lead  to  physical  understanding.  For 
example,  questions  of  how  a  structure  responds  to  various  excitations,  where  sources  of 


scattering  are  located,  or  how  discontinuities  affect  field  propagation,  can  be  answered  by 
simply  observing  electric  fields,  magnetic  fields,  voltages,  and/or  currents,  as  functions  of 
time  and  space. 

There  are  several  limitations  to  the  FDTD  method  and  these  limitations  are  often 
related  to  the  accurate  description  of  the  geometry  being  analyzed.  When  a  structure  has 
a  fine  feature  that  has  to  be  accurately  described,  a  fine  discretization  must  be  used  in  the 
computational  domain.  Complicated  geometries  are  also  difficult  to  model  because  all 
defining  dimensions  must  be  an  integral  multiple  of  the  cell  discretization  in  a  given 
direction.  This  can  lead  to  an  exorbitant  number  of  unknowns,  which  will  result  in  an 
intractable  problem  because  of  either  the  amount  of  memory,  or  the  amount  of  time, 
required  to  perform  the  numerical  simulation.  The  same  type  of  limitations  arise  when 
attempting  to  solve  large  problems.  Here,  large  refers  to  problems  on  the  order  of  18 
million  unknowns,  or  problems  with  dimensions  on  the  order  of  7  A,  by  7  A  by  7  A  where 
the  cell  discretization  is  on  the  order  of  A/20  and  the  platform  used  to  run  the  FDTD  code 
has  128  MB  of  RAM.  With  unlimited  computer  resources,  any  size  problem  can  be 
solved.  How  long  the  simulation  will  take  to  run,  however,  is  another  issue. 

The  problem  of  modeling  large  problems  or  geometries  with  fine  features  has 
been  addressed  by  several  researchers  who  have  reported  expansion  techniques  [1.4], 
sub-cell  gridding  [1.5],  and  sub-cell  modeling  [1.6].  Of  these  methods,  sub-cell  modeling 
is  the  most  accurate  and  can  be  incorporated  into  the  FDTD  update  equations  without  loss 
of  speed  and  does  not  require  temporal  interpolation.  Another  approach  for  solving  large 
problems  is  to  implement  the  FDTD  algorithm  on  a  parallel  architecture  [1.7].  This 
approach  does  permit  the  solution  of  large  problems,  but  can  be  expensive  in  terms  of 
simulation  times  because  of  the  need  to  communicate  large  amounts  of  data  at  every  time 
step. 

The  other  geometrical  feature  that  presents  difficulties  for  the  FDTD  method  is 
curves  or  angles  other  than  90°.  The  usual  method  for  modeling  these  types  of  features  is 


to  use  a  staircase  approximation  to  the  geometry.  Local  mesh  refinement  techniques  have 
also  been  reported  [1.8]  and  [1.9].  An  alternative  approach  is  to  use  the  curvilinear 
FDTD  method  [1.10].  This  method  does  a  good  job  at  accurately  modeling  geometries 
that  have  features  that  can  be  modeled  with  continuous  non -orthogonal  coordinate  lines. 
Special  conditions  and  update  equations  are  required  when  irregular  grids  are 
encountered.  The  accuracy  and  enhanced  modeling  capability  of  the  curvilinear  method 
certainly  does  not  come  without  a  price.  The  curvilinear  method  requires  at  least  twice  as 
much  memory  as  the  conventional  FDTD  method  to  store  covariant  and  contravariant 
field  components,  and  the  curvilinear  method  is  slower  because  there  are  more  operations 
in  the  update  equations  and  conversions  are  required.  In  many  instances,  both  the 
memory  requirements  and  computational  time  required  for  the  curvilinear  method  render 
the  method  prohibitively  expensive.  This  situation  should  be  remedied  to  an  extent  in  the 
future  by  faster  computer  processors  and  the  availability  of  machines  with  large  amounts 
of  memory. 

The  primary  goals  of  this  thesis  are  to  present  the  nonuniform  orthogonal  FDTD 
method,  and  to  demonstrate  its  applicability  to  a  variety  of  electromagnetic  problems. 
The  nonuniform  aspect  of  the  method  gives  great  flexibility  in  modeling  geometries. 
Furthermore,  the  variability  of  the  mesh  resolution  also  enables  the  user  to  move  the 
boundaries  of  the  computational  domain  farther  away  from  the  center  of  the  problem 
without  an  undue  increase  in  the  number  of  cells.  Most  significantly,  the  orthogonality  of 
the  method  preserves  the  speed  of  the  conventional  FDTD  method.  These  three  features 
of  the  nonuniform  orthogonal  FDTD  method  are  demonstrated  by  means  of  numerical 
examples  throughout  the  thesis. 

Long  simulation  run  times  can  also  be  an  undesirable  characteristic  of  the  FDTD 
method.  The  issue  of  temporal  truncation  and  extrapolation,  which  must  be  applied  with 
great  care,  is  not  addressed  in  this  thesis.  However,  numerical  techniques  which  enable 
the  extraction  of  frequency-dependent  information  from  a  minimum  number  of  FDTD 


simulations  are  presented.  Moreover,  these  same  techniques  can  be  applied  to  extract  the 
reflections  from  imperfect  absorbing  boundary  conditions. 

Since  the  simulation  time  is  directly  proportional  to  the  size  of  the  computational 
domain,  it  is  advantageous  to  keep  the  computational  domain  small.  This  involves  one  of 
the  most  important  issues  in  electromagnetic  modeling  using  the  finite  difference  time 
domain  method,  which  is  the  truncation  of  the  computational  domain  with  absorbing 
boundary  conditions  (ABCs).  An  underlying  theme  common  to  research  concerning 
ABCs  is  that  it  is  desirable  to  use  an  ABC  which  is  not  only  accurate  in  the  modeling  of 
waves  at  the  boundaries,  but  also  can  be  brought  close  to  the  radiating  source  or 
discontinuity,  thereby  saving  on  the  memory  and  computation  requirements.  Depending 
upon  the  problem  being  solved,  the  absorbing  boundary  should  simulate  the  outward 
propagation  of  traveling  waves  at  various  angles  of  incidence,  guided  modes,  or 
evanescent  modes.  Clearly,  it  is  advantageous  to  have  a  robust  absorbing  boundary 
condition  that  can  be  used  for  different  types  of  problems,  i.e.,  guided  wave  or  radiating 
wave.  The  modified  dispersive  absorbing  boundary  condition  is  such  an  ABC  [1.11],  and 
it  is  analyzed  and  implemented  on  a  nonuniform  grid  in  this  thesis. 

1.2  Outline  of  the  Thesis 

The  nonuniform  orthogonal  FDTD  method  is  applied  to  a  wide  variety  of 
electromagnetic  problems.  Whenever  possible,  the  FDTD  results  are  compared  with 
analytic  expressions,  experimentally  measured  results,  or  the  results  of  other  numerical 
methods,  including  results  from  circuit  simulators.  In  certain  instances,  due  to  the 
uniqueness  of  the  problem  analyzed,  comparisons  are  not  possible. 

Chapter  2  introduces  the  nonuniform  orthogonal  FDTD  method  after  briefly 
describing  the  conventional  FDTD  method.  The  nonuniform  algorithm  is  developed 
based  upon  the  general  curvilinear  FDTD  method.  Following  the  derivation  of  the 
nonuniform  orthogonal  update  equations,  the  error  due  to  the  nonuniform  grid  is 


discussed.  A  theoretical  analysis  is  followed  by  a  numerical  example.  Finally,  problems 
involving  the  finline  waveguide  and  lossy  layered  media  are  analyzed  using  the 
nonuniform  orthogonal  FDTD  method. 

Chapter  3  presents  the  implementation  of  the  dispersive  boundary  condition 
(DBC)  on  a  nonuniform  orthogonal  grid.  Motivation  for  the  study  of  the  DBC  is 
provided  with  an  analysis  of  the  Mur  ABCs.  The  DBC  and  modified  DBC  equations  are 
derived,  followed  by  a  discussion  on  stability.  The  DBC  for  nonuniform  grids  is  tested 
with  a  microstrip  line  problem  and  a  general  radiating  problem.  Error  analysis  is  carried 
out  in  the  time  domain,  the  frequency  domain,  and  spatially. 

Chapter  4  presents  a  variety  of  methods  for  calculating  the  frequency-dependent 
characteristics  of  microstrip  lines,  striplines,  and  discontinuities.  Use  of  a  Prony 
technique  with  the  FDTD  method  is  presented  and  applied  to  two-port  scattering 
parameter  calculations.  A  local  mesh  refinement  technique  is  introduced  for  triangular 
metallization  and  demonstrated  to  increase  the  accuracy  of  the  FDTD  method.  The 
FDTD  method  and  Prony  technique  are  combined  to  analyze  the  complicated  transition 
from  a  grounded  coplanar  waveguide  to  a  microstrip  line. 

Chapter  5  uses  the  FDTD  method  to  characterize  coupled  microstrip  lines  and 
coupled  striplines.  First,  the  symmetric  problem  is  considered  and  an  even  and  odd  mode 
theory  is  used  in  addition  to  the  methods  of  Chapter  4  to  characterize  the  lines.  The 
methods  of  Chapter  4  are  extended  to  handle  the  asymmetric  case.  This  new  method  is 
presented,  followed  by  numerical  examples  analyzing  coupled  symmetric  and  asymmetric 
lines.  A  multidielectric  asymmetric  coupled  stripline  problem  is  also  analyzed. 

Chapter  6  uses  the  FDTD  method  to  analyze  microstrip  antennas.  The  near-field 
to  far-field  transformation  is  discussed,  followed  by  the  analysis  of  a  complicated  balun 
fed  folded  dipole  antenna.  The  benefits  of  using  the  nonuniform  orthogonal  FDTD 
method  are  demonstrated  through  the  analysis  of  microstrip  patch  antennas.  The  results 
obtained  are  compared  with  experimental  measurements. 


Finally,  Chapter  7  summarizes  the  results  of  this  thesis  and  suggests  topics  for 


further  study. 
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CHAPTER  2 


THE  NONUNIFORM  ORTHOGONAL  FINITE  DIFFERENCE  TIME  DOMAIN 

METHOD 


2.1  Introduction 

In  this  chapter  the  Nonuniform  Orthogonal  Finite  Difference  Time  Domain 
method  is  presented.  After  first  considering  the  conventional  uniform  FDTD  method,  the 
nonuniform  algorithm  is  developed  based  upon  the  general  curvilinear  FDTD  method. 
The  development  entails  the  discussion  of  basis  vectors,  metrics,  and  field  locations. 
Following  the  derivation  of  the  nonuniform  orthogonal  update  equations,  the  error  due  to 
the  nonuniform  grid  is  discussed.  A  theoretical  analysis  is  followed  by  a  numerical 
example.  Finally,  examples  of  electromagnetics  problems  that  are  well-suited  to  the 
nonuniform  method  are  presented.  Specifically,  problems  involving  the  finline 
waveguide  and  lossy  layered  media  are  analyzed. 

2.2  Uniform  Finite  Difference  Time  Domain  Method 

The  Finite  Difference  Time  Domain  (FDTD)  method  was  introduced  to  the 
electromagnetics  community  by  Yee  in  1966  [2.1].  Because  of  its  versatility  and  the 
advances  in  computer  capability,  it  has  enjoyed  widespread  use  in  recent  years.  The 
FDTD  method  discretizes  Maxwell’s  equations  in  both  space  and  time  using  second- 
order  accurate  central  difference  formulas.  Arbitrary  geometries  are  described  on  a 
uniform  rectangular  mesh,  and  the  electric  and  magnetic  fields  are  determined  at  discrete 
locations  within  the  mesh  as  functions  of  time.  The  electric  field  values  are  located  on 
the  edges  of  the  rectangular  FDTD  cells,  and  the  magnetic  field  values  are  located  at  the 
centers  of  the  faces  of  the  cells.  The  dimensions  of  each  cell  are  Ax  by  Ay  by  Az,  and  the 
time  step  is  At.  The  time  step  is  related  to  the  mesh  discretization  through  the  Courant 


stability  criterion  shown  in  (2.1),  where  c  is  the  velocity  of  light  in  the  medium.  The 
method  becomes  unstable  when  (2.1)  is  violated,  but  for  optimal  performance,  At  should 
be  chosen  as  large  as  possible  [2.2]. 


A  t< — - - — ■  (2.1) 

n  i  r  (  J 

y  Ax2  Ay2  A z2 

The  electric  field  at  time  nAt  is  updated  in  terms  of  the  electric  field  at  the  same  location 
at  time  (n-l)At  and  the  surrounding  magnetic  fields  at  time  (n-l/2)At.  The  update 
equations  for  the  x-components  of  the  electric  and  magnetic  fields  are  given  in  (2.2)  and 
(2.3). 
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where  Fn(i,j,k)  =  F(iAx,jAy,kAz\nAt),  c  is  the  electrical  conductivity,  |i  is  the 
permeability,  and  e  is  the  permittivity.  Update  equations  for  the  other  field  components 
are  obtained  by  permuting  the  fields  and  indices. 


2.3  Nonuniform  Orthogonal  Finite  Difference  Time  Domain  Method 

In  this  section,  the  update  equations  for  nonuniform  orthogonal  FDTD  will  be 
derived  from  curvilinear  FDTD  update  equations.  Before  deriving  the  update  equations, 
the  basis  vectors,  metrics,  covariant  fields,  contravariant  fields,  and  physical  fields  will  be 
discussed.  The  presentation  here  will  follow  that  of  Holland  [2.3].  Any  general 
coordinate  system  can  be  characterized  by  unitary  basis  vectors,  (a,,^,^).  In  general, 
these  vectors  need  not  be  orthogonal  to  each  other  and  are  not  of  unit  length.  However, 
for  the  purposes  of  this  work,  the  vectors  will  be  orthogonal  to  each  other.  Typical  basis 
vectors  used  in  this  work,  shown  in  Fig.  2.1,  are  parallel  to  the  edges  of  the  cell. 
Reciprocal  basis  vectors  can  also  be  used  to  characterize  the  coordinate  system  and  can 
be  defined  in  terms  of  the  basis  vectors  as  follows: 

a'  =  aj  xak/^[g  (2.4) 

where  Jg  is  equal  to  the  volume  of  the  rectangular  prism  formed  by  the  basis  vectors 
(a,,<32,a3).  These  orthogonal  vectors  satisfy  the  following  properties: 

(2.5) 


Figure  2.1  Basis  vectors  shown  at  two  points  on  a  typical  FDTD  cell.  Point  1  is  located 
at  the  center  of  an  edge  along  the  u 1  coordinate  line,  and  Point  2  is  located 
at  the  center  of  a  face  normal  to  the  u 1  coordinate  line. 
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3,Sj-8tg'  (2.6) 

(2.7) 

where  gy  is  the  inverse  of  g,J  and  6y  is  the  Kronecker  delta  function. 

A  general  vector  can  be  written  in  terms  of  its  covariant  components,  as  in  (2.8), 
or  in  terms  of  its  contravariant  components,  as  in  (2.9). 
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i=l 

(2.8) 

i=l 

(2.9) 

ng  formulas: 

e  =  p..ej8 

l  OIJ  IJ 

(2.10) 

e‘=giJejSy 

(2.11) 

The  e{  and  e'  are  not  physical  fields.  The  physical  fields  are  related  to  the  et  and 
e'as  follows: 

Ei=4¥ei  (2.12) 

and  the  physical  fields  are  related  by  the  following  expressions 


where 


E  =  G  EJ 

l  ij 

E‘  =  GijE: 


Gij  =  -&!rgiJ 
\gJJ 


(2.14) 

(2.15) 

(2.16) 

(2.17) 


However,  from  the  above  definitions  for  g0  and  g‘J ,  (2.14)  and  (2.15)  simplify  to 

Using  the  above  definitions,  the  first  component  of  Faraday’s  law  in  point  form 


can  be  written  as 
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(2.18) 


where  the  u‘  are  the  coordinate  lines  of  the  curvilinear  system  and  the  permeability  is  not 
a  function  of  time.  Similarly,  the  first  component  of  Ampere’s  law  in  point  form  can  be 
written  as 


dE1  ]  (77  H 

eir+aE  Jf  a? 


s[hJ4?)  s[hJ4?) 

du2  du 3 


(2.19) 


A  leapfrog  scheme  is  usually  used  to  solve  (2.18)  and  (2.19).  First,  (2.18)  is  solved  for 
the  H\  which  are  then  converted  to  Hr  Then,  (2.19)  is  solved  for  the  £".  The  £"  are 
converted  to  Ei  and  used  to  solve  (2.18),  and  the  process  is  repeated.  For  the  orthogonal 
case,  it  has  been  shown  that  Et  =  E‘ .  Consequently,  the  conversions  are  not  necessary. 
Moreover,  it  will  be  demonstrated  that  the  discretized  versions  of  (2.18)  and  (2.19)  reduce 
to  simple  forms  requiring  the  same  number  of  operations  as  (2.2)  and  (2.3). 

Using  central  differences,  (2.18)  is  written  as 

H'(ij,ky+V2 = H\i,j,k)n~y2 +Si- 


f  E2(i,j,k)‘  /4fF  -  E2(i,j,k  -l)’  /4?* 
y  U.{t,j,k)-U.(i,j,k- 1) 

E3(iJ-l,k)-/4^-E,(i,j,kT/4g» 
U.(i,i,k)-U,(i,j- U) 


(2.20) 


where  fi  and  ^gu/g  are  evaluated  at  the  Hx  mesh  points  and  the  ^[g^  are  evaluated  at 
the  Et  mesh  points,  as  the  notation  indicates.  The  Un  s  are  the  successive  points  of  the 
curvilinear  coordinates  and,  by  definition,  their  difference  is  unity. 


Since  the  nonuniform  orthogonal  system  considered  in  this  work  is  Cartesian,  the 
coordinate  lines  (w',m2,m3)  correspond  directly  to  the  x-,  y-,  and  z-directions.  Figure  2.2 
shows  that  the  Ex  and  Hx  ( El  and  Hx )  components  are  located  at  positions  marked  with 
unfilled  and  filled  triangles,  respectively.  The  y-  and  z-component  positions  are  marked 
with  circles  and  squares,  respectively.  Cell  nodes  are  located  at  cell  vertices  and  labeled 
with  Un ,  and  cell  face  centers  are  labeled  with  Uc.  The  following  definitions  are  made  to 
assist  in  the  simplification  process.  Let  dxn(i)  be  the  distance  between  two  nodes  located 
at  i  and  i-1.  Similarly,  let  dyn(j)  be  the  distance  between  two  nodes  at  j  and  j-1,  and 
dzn(k )  be  the  distance  between  nodes  at  k  and  k-1.  Define  dxc(i)  as  the  distance 
between  the  cell  centers  U2c{i+\,j,k)  and  U2{i,j,k)  as  shown  in  Fig.  2.2.  Note  that 
dxc(i)  is  also  equal  to  the  distance  between  cell  centers  U3c(i+l,j,k)  and  U3c(i,j,k).  Let 
dyc(j)  and  dzc(k)  be  defined  similarly. 


Figure  2.2  Typical  FDTD  cell  showing  locations  of  the  electric  and  magnetic  fields. 

Electric  fields  are  located  along  cell  edges  and  magnetic  fields  are  located  at 
cell  face  centers,  labeled  Ur.  The  cell  nodes  are  labeled  U  . 

c  n 
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Using  these  definitions,  the  various  factors  in  (2.20)  can  be  evaluated  and 


simplified. 


Sn\H]{iJ,k)=dxc(i)-dxc(i) 

= ■  dyn(j)  ■ dzn^ 

5=1  *  V? I  Vdyn{j) 

^2(i,J,k)  E2(i,j,k- 1) 


V?|  “V?  =  1/  dzn(k) 


(2.21) 


(2.22) 


(2.23) 


(2.24) 


Substituting  (2.21)  through  (2.24)  into  (2.20)  and  simplifying  result  in  the  following 
equation. 

H'{iJ,k)n+v 2  =  H'{i,j,k)n-y2+  — 


E2  O',  j\  kf  -  E2  (i,  j,  k-iy  |  £3  (i,  l  - 1,  fc)"  -  £3  (i,  j,  k)n 


dzn(k) 


dyn{j ) 


Equation  (2.25)  is  identical  in  form  to  (2.3). 

Turning  now  to  (2.19),  central  differencing  leads  to  (2.26). 


(2.25) 


— —  f f  At)  /(  oAt 

Sl/*  7/1+¥ 


#2  (iJ,k)H*V2/ ~  H2(i,j,k  +  l)n+1/2/ Vp1 
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#3(*  J  + 1  ,A:)B+V2/ VP1  -  H2(iJ,k)n+l/2/4^ 
U3c{i,j  +  l,k)-U3c{i,j,k) 


(2.26) 


Again,  note  that  the  denominators  of  the  partials  with  respect  to  the  curvilinear 
coordinates  are  unity  by  definition.  The  metric  factors  are  evaluated  in  the  following 
equations. 
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8n\El{u,k)=dxn(i)-dxn(i) 

Me ,(i,m  =  ■  dyc(j) '  dzc(^> 


S'" I  f...  =vr*  =  l/dyc(j) 

#2  (*»,/»*)  ;,*+!) 


V?  =  V?1  =  l/dzc(k) 

Substituting  Equations  (2.27)  through  (2.30)  into  (2.26)  results  in 

(Y,  crAr^  tf.  .  crAtMpw.  .  .y, 
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(2.27) 

(2.28) 

(2.29) 

(2.30) 


+ 
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7/3(i,7  +  l,fc)”+1/2-//3(t-,y,fc) 
dyc{j ) 


«+l/2  \\ 


77 


(2.31) 


which  is  identical  in  form  to  Equation  (2.2). 

It  is  important  to  note  that  since  the  nonuniform  orthgonal  update  equations  are 
identical  in  form  to  the  uniform  update  equations,  there  is  no  loss  of  speed  when  updating 
the  fields.  Care  must  be  taken  when  computing  the  coefficients  in  the  nonuniform  update 
equations  since  variable  cell  discretization  must  be  taken  into  account  in  order  to 
accurately  weight  the  coefficients  when  dealing  with  interfaces  between  different  media 
[2.4].  However,  the  advantages  from  memory  savings  and  decrease  in  total  run  time  due 
to  a  smaller  number  of  unknowns  far  outweigh  this  slight  disadvantage.  For  many  classes 
of  problems,  this  difficulty  can  be  avoided  entirely  by  requiring  that  cell  discretizations 
normal  to  a  material  interface  be  equal,  in  which  case,  the  algorithm  can  be  implemented 
with  the  same  efficiency  as  the  uniform  method. 
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2.4  Error  from  Grid  Dispersion  on  Nonuniform  Grids 

It  is  clear  that  the  ability  to  allow  the  FDTD  mesh  to  vary  will  result  in  a  reduction 
in  the  number  of  unknowns.  The  mesh  can  be  dense  in  areas  where  fine  features  must  be 
accurately  modeled  and  then  allowed  to  expand  or  grow  to  a  relatively  "coarse"  mesh  in 
areas  where  fine  discretization  is  not  necessary.  However,  one  cannot  simply  allow  the 
mesh  to  grow  at  an  arbitrary  rate.  The  mesh  must  be  varied  gradually;  otherwise,  phase 
error  due  to  grid  dispersion  will  contaminate  the  solution  at  an  unacceptable  level.  An 
excellent  discussion  of  grid  dispersion  error  can  be  found  in  [2.5].  The  development  in 
[2.5]  will  be  followed  here,  simplified  to  the  one-dimensional  case  and  perturbed  to 
include  the  nonuniformity  of  the  grid. 

Consider  the  scalar  wave  equation  for  homogenous  media  given  in  (2.32). 

^-j^(z,O  =  V20M  (2.32) 

Since  any  wave  can  be  expressed  as  a  sum  of  plane  waves,  solutions  of  (2.32)  are 
comprised  of  plane  waves. 

(p(z,t)  =  A(t)e-Jk‘z  (2.33) 

The  discretized  form  of  the  plane  wave  in  (2.33)  is  given  in  (2.34). 

</>‘m  =  **%“*"*  (2.34) 

Discretizing  (2.32)  using  central  differencing  results  in  (2.35). 

€'  =  Jc2  +<-,)  +  2*L  -  '  (2.35) 

Substituting  (2.34)  into  (2.35)  results  in  the  one-dimensional  numerical  dispersion 
relation  given  in  (2.36). 


sin 


(  (oAt 


cAt\  .  ( k.Az\ 

AzJ  l  2 


(236) 


Now,  consider  a  cell  whose  discretization,  A z2,  is  larger  than  A z-  Solving  (2.36) 


for  the  numerical  propagation  constant  yields 
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j3  =  —  sin"1 
At 


lAz> 


(2.37) 


where  f5  =  (0/c  and  k0-kz.  Equation  (2.37)  is  an  exact  solution  to  the  numerical 
dispersion  relation  (2.36).  An  approximate  solution  to  (2.36)  can  be  obtained  by  writing 


.  ( a>At 
sin  - 

l  2 


=  rs 


(2.38) 


where  r  =  ~^  a°d  s  =  sin^^^j.  Using  series  expansions  for  sin*  and  sin"1  x  and 
retaining  second-order  terms  yield  the  approximate  solution  to  the  numerical  dispersion 
relation 


P  =  K 


{  (K^i)  (cAtfk, 

24  24 


(2.39) 


where  f3  =  co/c  and  k0  =  kz. 

Suppose  that  A z,  is  the  minimum  cell  discretization  and  that  At  is  chosen  such 
that  At  =  AzJc.  Clearly,  there  is  no  phase  error  due  to  grid  dispersion  if  A z2  =  Az,  in 
(2.37)  and  (2.39).  When  A z2  *  A z, ,  then  (2.37)  and  (2.39)  can  be  used  to  determine  the 
normalized  phase  error  for  a  given  cell  as  a  function  of  the  ratio  Az2  /  Az, .  For  example, 
let  Az,  =  10  mm  and  allow  A z2  to  vary.  At  a  frequency  of  1.5  GHz,  Az,  =  A/20.  Since 
the  largest  practical  value  of  Az2  is  A  / 10,  the  ratio  Az2  /  Az,  can  range  from  1  to  2.  The 
percent  error  in  normalized  phase  is  plotted  in  Fig.  2.3,  where  Ml  refers  to  (2.37)  and  M2 
refers  to  (2.39)  and  the  frequency  values,  f,  are  given  in  gigahertz.  The  agreement 
between  methods  Ml  and  M2  verifies  the  approximations  made  in  deriving  (2.39).  For 
f  =  1 .5  GHz,  there  is  a  1 .2%  error  in  the  normalized  phase  at  a  ratio  of  2.  For 
f  =  0.75  GHz,  the  ratio  goes  up  to  4  where  the  error  is  about  1.55%.  This  error  is  per  cell, 
so  the  error  will  accumulate  as  the  field  propagates  from  cell  to  cell.  It  is  evident  that  the 
smaller  the  cell  size  in  terms  of  wavelength,  the  larger  the  acceptable  growth  ratio 
between  cells.  For  example,  there  is  less  than  0.2%  error  at  f  =  1.5  GHz  for  a  ratio  of  1.2 
and  less  than  0.2%  error  at  f  =  0.75  GHz  for  a  ratio  of  1.6. 
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Figure  2.3  Percent  error  in  normalized  phase  per  cell  for  two  frequencies  and  two 
methods.  Ml  refers  to  (2.37)  and  M2  refers  to  (2.39). 


In  practical  applications,  At  is  usually  chosen  less  than  the  value  specified  by  the 
Courant  stability  criterion,  (2.1).  To  show  how  this  affects  the  error  due  to  the 
nonuniform  grid,  consider  A z,  =  10  mm  and  f  =  1.5  GHz.  Rewriting  (2.39)  leads  to 


P  =  K 


(2.40) 


where  =  co/c,  ka=kz,  and  cAt  =  r\Azx.  Figure  2.4  shows  that  varying  the  time  step 
increases  the  percent  error  in  the  normalized  phase  by  not  more  than  0.1  %  when  T]  is  at 
least  0.9. 

To  test  the  effect  of  the  nonuniformity  of  the  grid  on  a  practical  numerical 
example,  a  stripline  embedded  in  free  space  was  analyzed.  Since  the  stripline  structure 
supports  a  TEM  mode,  the  grid  dispersion  error  can  be  determined  by  varying  the  grid  in 
the  direction  of  propagation  and  calculating  the  relative  dielectric  constant.  In  this  case, 
the  relative  dielectric  constant  should  be  one.  The  geometry  of  the  stripline  is  shown  in 
Fig.  2.5(a),  where  the  separation  between  the  ground  planes,  H,  is  6.0  mm,  and  the  width, 


Figure  2.4  Percent  error  in  normalized  phase  per  cell  for  f  =  1.5  GHz  and  77  =  cAtl  Az, 
where  Az,  =  A/20. 


W,  of  the  centered  stripline  is  4.0  mm.  The  cross-section  of  the  stripline  is  discretized 
uniformly  with  Ax  =  Ay  =  500  |im.  The  line  runs  in  the  z-direction,  and  an  example  of 
15  cells  of  uniform  cell  discretization  is  shown  in  Fig.  2.5(b),  and  a  nonuniform 
discretization  with  growth  rate  1.2  is  shown  in  Fig.  2.5(c).  Note  that  in  Fig.  2.5(c)  the 
mesh  begins  expanding  with  the  12th  cell  and  stops  at  the  15th  cell.  The  first  eleven  cells 
are  uniformly  discretized.  The  smallest  cell  in  the  z-direction  is  Az]  =  500  (im  and  the 
largest  cell  is  A z2  =  1.0  mm. 

The  effective  dielectric  constant  is  calculated  by  monitoring  the  current  at  two 
locations  along  the  line,  15  cells  apart.  Depending  on  the  cell  discretization,  the  distance, 
L,  between  the  two  points  will  vary.  The  time  domain  signatures  of  the  currents  at  the 
two  positions  are  transformed  into  the  frequency  domain  and  the  effective  dielectric 
constant  is  then  given  by  (2.41) 


In 


'  M  v 

I(d  +  L) 


-l2 


(2.41) 


where  /  is  a  function  of  frequency. 
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Figure  2.5  Geometry  of  stripline  used  to  test  numerical  grid  dispersion,  (a)  Cross- 
section  of  stripline,  W  =  4.0  mm,  H  =  6.0  mm.  (b)  Uniform  cell 
discretization  in  the  direction  of  propagation,  (c)  Nonuniform  cell 
discretization  in  the  direction  of  propagation  with  growth  rate  1.2,  cells 
growing  near  the  "d+L"  location. 

Six  different  growth  rates  were  tested  and  the  value  of  the  effective  dielectric 
constant  is  shown  in  Fig.  2.6,  where  the  curves  are  labeled  on  the  graph  according  to  the 
growth  rate  of  the  mesh.  Note  that  the  uniform  case  (growth  rate  =  1.0)  results  in  a 
maximum  error  of  0.75  %  at  30  GHz.  The  nonuniform  cells  were  positioned  within  the 
15-cell  region  such  that  the  largest  cell  was  the  15th  cell.  This  explains  the  fact  that  a 
growth  rate  of  1.1  gives  rise  to  more  total  error  than  larger  growth  rates  for  frequencies 
up  to  about  23  GHz.  When  the  growth  rate  is  1.1,  it  takes  seven  cells  to  go  from  Az,  to 
Az2,  whereas  using  larger  growth  rates,  a  much  smaller  number  of  cells  is  used  in  the 
transition  from  Az,  to  Az2.  When  the  mesh  grows  from  the  first  monitoring  location,  the 
total  error  exhibits  the  expected  behavior  as  a  function  of  frequency  and  the  growth  rate, 


as  shown  in  Fig.  2.7.  Error  accumulates  as  the  field  propagates  from  "d"  to  "d+L"  and  the 
larger  cells  contribute  more  error  per  cell.  For  a  growth  rate  of  1.2,  the  accumulated  error 
is  about  1.32  %  at  30  GHz. 


Figure  2.6  Effective  dielectric  constant  as  a  function  of  frequency  for  the  air-filled 
stripline  of  Fig.  2.5.  The  growth  rate  varies  as  indicated  in  the  legend,  and 
the  mesh  grows  near  the  second  monitoring  point,  "d+L." 
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Figure  2.7  Effective  dielectric  constant  as  a  function  of  frequency  for  the  air-filled 
stripline  of  Fig.  2.5.  The  growth  rate  varies  as  indicated  in  the  legend,  and 
the  mesh  grows  near  the  first  monitoring  point,  "d." 

2.5  Finline  Waveguide 

The  finline  waveguide  became  popular  in  the  late  1970s  with  the  increased 
interest  in  millimeter- wave  integrated  circuits  [2.6].  The  finline  structure  is  shown  in 
Fig.  2.8.  It  can  be  viewed  as  a  shielded  slotline,  a  ridged  waveguide  with  dielectric,  or  a 
slab  loaded  waveguide  with  fins  [2.7].  For  millimeter-wave  applications,  the  finline  is 
preferred  over  the  microstrip  line  because  the  finline  does  not  require  stringent 
manufacturing  tolerances,  is  less  susceptible  to  the  propagation  of  higher-order  modes, 
and  interfaces  well  with  waveguide  instrumentation  [2.6].  Several  frequency  domain 
techniques  have  been  succesfully  applied  to  characterize  the  finline  waveguide  and 
related  structures  [2.7]-[2.11].  This  type  of  problem  is  particularly  well-suited  to 
nonuniform  orthogonal  FDTD  because  of  the  thin  dielectric  slab  which  has  to  be  modeled 
accurately.  The  mesh  can  grow  as  it  moves  away  from  the  gap  and  toward  the  shield 


walls. 
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Figure  2.8  Geometry  of  the  unilateral  finline  waveguide  (with  dielectric)  and  the 
symmetrical  ridged  waveguide  (without  dielectric). 


In  this  section,  the  cutoff  frequency  of  a  symmetrical  ridged  waveguide  with  a 
varying  aperture  is  calculated.  The  geometry  being  investigated  is  shown  in  Fig.  2.8, 
where  the  dielectric  has  a  dielectric  constant  of  one.  The  entire  cross-section  is  excited 
with  an  impulse  source,  and  the  time  signature  of  Ey  is  monitored  at  several  locations. 
These  time  signatures  are  then  Fourier  transformed  to  the  frequency  domain,  and  a  spike 
occurs  at  the  cutoff  frequency.  The  aperture  width,  w,  was  varied  and  is  shown  as  a 
fraction  of  b,  the  length  of  the  top  wall  of  the  waveguide,  as  shown  in  Fig.  2.8. 
Normalized  cutoff  frequencies  calculated  with  the  FDTD  method  are  compared  with 
results  using  the  regular  solution  of  the  singular  integral  equation  (RSSIE)  method  [2.11]. 
As  shown  in  Fig.  2.9,  the  agreement  between  the  two  methods  is  very  good. 

After  monitoring  fields  in  the  time  domain  and  transforming  to  the  frequency 
domain,  the  field  distribution  is  readily  available  over  a  wide  frequency  band.  Figure 
2.10  shows  the  normalized  distribution  of  Ey  in  the  aperture  of  a  unilateral  finline 
waveguide  with  dielectric  constant  2.22  at  a  frequency  of  30  GHz  for  two  different 
aperture  widths.  Again,  the  entire  cross-section  is  excited  with  an  impulse  and  the  time 


signatures  are  monitored  in  the  aperture  and  then  transformed  to  the  frequency  domain. 
These  results  compare  well  with  aperture  distributions  calculated  using  the  spectral 
domain  method. 


□  FDTD  ■  RSSIE 


1  0.6667  0.5  0.3333  0.25  0.125 
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Figure  2.9  Normalized  cutoff  frequency  of  symmetrical  ridged  waveguide  computed 
using  FDTD  and  RSSIE. 


Figure  2.10  Distribution  of  normalized  Ey  in  aperture  of  unilateral  finline  waveguide  for 
two  different  aperture  widths. 
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2.6  Lossy  Layered  Media 

An  application  that  is  ideal  for  nonuniform  FDTD  is  the  analysis  of  fields  in  a 
lossy  medium.  The  lossy  medium  considered  in  this  section  is  seawater,  which  is 
assumed  to  have  a  constant  conductivity  of  4  S/m  and  er  =  72.  In  sediment  layer  #1, 
cr,  =  2  S/m  and  in  sediment  layer  #2,  a2  =  1  S/m.  The  geometry,  shown  in  Fig.  2.1 1,  is  a 
10  m  long  dipole  oriented  along  the  x-axis  and  excited  with  a  Gaussian  pulse  that  has  a 
3  dB  cutoff  frequency  of  10  kHz.  The  dipole  is  centered  in  a  volume  that  is  200  m  by 
100  m  by  100  m.  The  problem  is  discretized  with  a  mesh  that  has  Ax  =  2  m,  Ay  =  1  m, 
and  Az  =  1  m  near  the  dipole.  The  mesh  discretization  then  expands  away  from  the 
source,  and  is  truncated  with  pec  boundaries.  The  time  step  for  FDTD  simulations  is 
determined  by  the  Courant  stability  criterion  as  shown  in  (2.1).  For  this  problem,  the 
frequencies  of  interest  are  extremely  low;  thus,  the  time  step  would  be  restrictively  small. 
Therefore,  a  different  approach  is  taken  for  describing  the  medium  and  choosing  the  time 
step. 


Figure  2.11  Geometry  used  for  the  analysis  of  fields  in  lossy  layered  media.  The  line  of 
observation  points  is  14  m  below  the  dipole  source.  The  dipole  source  is 
centered  in  the  computational  volume. 
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Since  c  =  \/^noiLr£0£r  ,  the  time  step,  At,  will  increase  if  er  increases.  The 
inequality  (oe«a  must  also  be  satisfied  because  seawater  is  a  good  conductor.  Since 
the  frequencies  of  interest  are  very  low,  er  can  be  increased  to  a  value  on  the  order  of 
1.0- 106  without  affecting  the  low  frequency  results  [2.12].  This  results  in  a  much  larger 
time  step.  The  value  of  er  is  chosen  to  be  4.0  •  106  and  At  is  chosen  to  be  2.2  |xs.  Electric 
and  magnetic  fields  are  monitored  along  a  line  parallel  to  the  x-axis  14  m  below  the 
dipole  source.  The  simulation  is  run  for  6000  iterations,  with  field  values  recorded  every 
ten  iterations.  The  time  signatures  are  converted  to  the  frequency  domain  via  the  discrete 
Fourier  transform  (DFT). 

When  the  entire  background  medium  is  seawater,  an  analytic  expression  for  the 
fields  exists.  Nonuniform  FDTD  results  are  compared  to  analytic  expressions  in 
Figs.  2.12  -  2.14.  The  agreement  is  excellent. 


Figure  2.12  Magnitude  of  Ex  in  seawater  observed  14  m  below  the  source  at  f  =  1  Hz. 


Position  on  Observation  Line  (m) 


Figure  2.13  Magnitude  of  Ez  in  seawater  observed  14  m  below  the  source  at  f  =  1  Hz. 


Figure  2.14  Magnitude  of  By  in  seawater  observed  14  m  below  the  source  at  f  =  1  Hz. 
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When  layers  of  sediment  are  present,  as  shown  in  Fig.  2.11,  the  FDTD  results  are 
compared  to  those  for  an  integral  equation  method  provided  by  the  U.S.  Navy  [2.13]. 
Figures  2.15  -  2.17  show  the  nonuniform  FDTD  and  Navy  results,  and  again  the 
agreement  is  excellent. 


Position  on  Observation  Line  (m) 

Figure  2.15  Magnitude  of  Ex  in  seawater  with  lossy  layers  observed  14  m  below  the 
source  at  f  =  1  Hz. 


Position  on  Observation  Line  (m) 


Figure  2.16  Magnitude  of  Ez  in  seawater  with  lossy  layers  observed  14  m  below  the 
source  at  f  =  1  Hz. 
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Position  on  Observation  Line  (m) 

Figure  2.17  Magnitude  of  By  in  seawater  with  lossy  layers  observed  14  m  below  the 
source  at  f  =  1  Hz. 

Figure  2.18  shows  the  magnitude  of  the  residual  electric  field,  which  is  the 
absolute  value  of  the  difference  between  the  fields  in  seawater  alone  and  the  fields  with 
the  layers  of  sediment  present.  The  agreement  between  the  results  of  the  two  methods  is 
excellent. 

The  computational  volume  used  for  the  FDTD  simulations  was  62  cells  by  62 
cells  by  62  cells.  When  the  same  problem  was  solved  using  uniform  FDTD,  the 
computational  volume  was  100  cells  by  100  cells  by  100  cells.  The  nonuniform  FDTD 
resulted  in  a  76%  reduction  of  required  memory.  Moreover,  the  smaller  number  of 
unknowns  also  led  to  a  shorter  computation  time.  To  further  demonstrate  the  utility  of 
the  nonuniform  mesh,  the  seawater  problem  was  solved  with  a  computational  volume  that 
was  39  cells  by  39  cells  by  39  cells.  Field  magnitudes  are  plotted  in  Figs.  2.19  -  2.21. 
Again  the  agreement  is  excellent.  This  FDTD  simulation  was  nine  times  faster  than  the 


62  by  62  by  62  case  and  gave  a  75%  memory  savings.  Relative  to  the  100  by  100  by  100 
case,  the  39  by  39  by  39  case  represents  a  94%  savings  in  memory. 


Ex  R  Navy  O  -  Ex  Nonuniform  FDTD 

Ez  R  Navy  Ez  Nonuniform  FDTD 


Position  on  Observation  Line  (m) 


Figure  2.18  Comparison  of  FDTD  results  with  integral  equation  results  for  the  residual 
electric  field  observed  14  m  below  the  source  at  f  =  1  Hz. 


Figure  2.19  Magnitude  of  Ex  in  seawater  observed  15.9  m  below  the  source  at  f  =  1  Hz. 


Position  on  Observation  Line  (m) 

Figure  2.20  Magnitude  of  Ez  in  seawater  observed  14.8  m  below  the  source  at  f  =  1  Hz. 


Position  on  Observation  Line  (m) 


Figure  2.21  Magnitude  of  Hy  in  seawater  observed  14.8  m  below  the  source  at  f  =  1  Hz. 
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The  nonuniform  FDTD  method  is  very  applicable  to  this  type  of  problem  because 
of  the  ability  to  let  the  mesh  discretization  grow.  Moreover,  a  single  FDTD  simulation 
provides  a  wide  band  of  frequency  information,  provided  the  geometry  is  excited  with  a 
transient  source.  This  was  the  case  here,  although  results  are  only  shown  for  a  single 
frequency.  The  nonuniform  orthgonal  FDTD  approach  is  preferred  over  moment  method 
type  solutions  for  this  problem  because  of  the  ability  to  vary  the  geometrical  and 
electrical  features  of  the  problem  without  the  derivation  of  a  complicated  Green's 
function.  For  example,  it  would  be  a  simple  matter  to  vary  the  surface  of  the  sediment  in 
a  random  fashion  so  as  to  model  a  rough  surface,  or  to  embed  any  kind  of  dielectric  or 
metal  slabs  horizontally,  or  vertically,  within  the  layers  of  sediment,  or  in  the  seawater. 

2.7  Conclusions 

This  chapter  described  the  uniform  FDTD  method  and  introduced  the  nonuniform 
orthogonal  FDTD  method.  The  nonuniform  orthogonal  FDTD  method  was  developed 
from  the  general  curvilinear  FDTD  method.  The  grid  dispersion  error  due  to  the 
nonuniform  mesh  was  analyzed,  theoretically  in  terms  of  phase  error  per  cell,  and 
numerically  in  terms  of  accumulated  phase  error  from  the  field  propagating  along  a 
nonuniformly  meshed  stripline.  It  was  demonstrated  that  the  cell  size  in  terms  of 
wavelength  as  well  as  the  growth  rate  of  the  mesh  are  critical  in  determining  the  amount 
of  error,  which  can  be  maintained  at  acceptably  low  levels  while  still  reducing  memory 
requirements.  Several  example  problems  were  presented  which  demonstrated  the 
flexibility  of  the  nonuniform  FDTD  method  as  well  as  the  attractive  features  of  memory 
savings  and  accurate  fine  feature  modeling.  In  the  case  of  the  finline  problem,  small  cell 
discretization  was  used  to  model  the  fine  features  of  the  geometry.  In  the  case  of  the 
lossy  layered  media  problem,  the  mesh  grew  rapidly  from  the  regions  of  interest,  which 
enabled  the  boundaries  of  the  computational  domain  to  be  positioned  sufficiently  far 
away  from  the  source  and  layers  of  sediment.  This  was  accomplished  with  a  small 
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number  of  unknowns  relative  to  the  uniform  method.  The  results  presented  were  in 
excellent  agreement  with  analytic  and  moment  method  techniques. 
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CHAPTER  3 

THE  DISPERSIVE  ABSORBING  BOUNDARY  CONDITION  APPLIED  TO 
NONUNIFORM  ORTHOGONAL  MESHES 


3.1  Introduction 

One  of  the  most  important  issues  in  electromagnetic  modeling  using  the  finite 
difference  time  domain  method  is  the  truncation  of  the  computational  domain  with 
absorbing  boundaries.  Considerable  efforts  have  been  directed  to  the  development, 
analysis,  and  improvement  of  absorbing  boundary  conditions  (ABCs)  [3. l]-[3. 15].  An 
underlying  theme  common  to  research  concerning  ABCs  is  that  it  is  desirable  to  use  an 
ABC  which  is  not  only  accurate  in  the  modeling  of  waves  at  the  boundaries,  but  also  can 
be  brought  close  to  the  radiating  source  or  discontinuity,  thereby  saving  on  the  memory 
and  computation  requirements.  Depending  upon  the  problem  being  solved,  the  absorbing 
boundary  should  simulate  the  outward  propagation  of  travelling  waves  at  various  angles 
of  incidence,  guided  modes,  or  evanescent  modes.  Clearly,  it  is  advantageous  to  have  a 
robust  absorbing  boundary  condition  that  can  be  used  for  different  types  of  problems,  i.e., 
guided  wave  or  radiating  wave.  The  modified  dispersive  absorbing  boundary  condition  is 
such  an  ABC  [3.6]. 

The  purpose  of  this  chapter  is  to  present  the  implementation  of  the  dispersive 
boundary  condition  (DBC)  on  a  nonuniform  orthogonal  grid.  Motivation  for  the  study  of 
the  DBC  is  presented  with  an  analysis  of  the  Mur  ABCs.  The  DBC  and  modified  DBC 
equations  are  derived,  followed  by  a  discussion  on  stability.  The  DBC  for  uniform  grids 
is  numerically  tested  with  a  microstrip  line  problem.  Finally,  the  DBC  for  nonuniform 
grids  is  tested  with  a  microstrip  line  problem  and  a  general  radiating  problem.  Error 
analysis  is  carried  out  in  the  time  domain,  the  frequency  domain,  and  spatially. 


3.2  The  Mur  Absorbing  Boundary  Conditions 

It  is  worthwhile  to  examine  the  limitations  of  the  first-  and  second-order  Mur 
absorbing  boundary  conditions  [3.1].  Not  only  does  this  motivate  the  study  of  the 
dispersive  boundary  condition  (DBC),  but  it  also  provides  some  insight  regarding  the 
performance  of  the  DBC. 

The  well-known  first-order  Mur  ABC  is  given  in  Equation  (3.1). 


This  ABC  is  unconditionally  stable,  straightforward  to  implement,  and  has  minimal 
memory  requirements.  It  can  be  used  for  radiating  problems  as  well  as  for  guided  wave 
problems.  The  main  drawback  of  the  first-order  Mur  is  that  it  only  absorbs  waves  that  are 
normally  incident  with  velocity  c.  Additionally,  higher-order  ABCs  are  usually  preferred 
over  the  first-order  Mur  because  the  reflections  from  the  first-order  Mur,  which  are 
typically  on  the  order  of  3  to  4%  of  the  field  incident  on  the  boundary,  can  seriously 
degrade  frequency  domain  results. 

The  second-order  Mur  ABC  is  given  in  Equation  (3.2). 


Id2  Id2  if  d2  d2 
cdxdt  c2  dt2+ 2{dy2+ dz2 


ELo=0 


(3.2) 


The  second-order  Mur  requires  tangential  derivatives  on  the  boundary.  Obviously,  this 
leads  to  problems  on  the  comers  of  the  computational  domain  since  there  is  not  enough 
information  to  perform  the  derivative.  Thus,  additional  comer  conditions  must  be 
implemented.  A  more  serious  concern  is  the  inability  of  the  second-order  Mur  to  absorb 
guided  waves,  for  example,  waves  traveling  on  a  microstrip  line.  This  phenomenon  can 
be  explained  by  considering  the  dispersion  relation  and  the  Mur  operator. 

The  dispersion  relation  is  given  in  Equation  (3.3) 
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where 


kx  =  fcsin(0)cos(0) 

(3.4) 

ky  =  fcsin(0)sin(0) 

(3.5) 

kz  =  kcos(6) 

(3.6) 

are  the  Cartesian  components  of  the  propagation  vector,  6  is  measured  from  the  z-axis 
and  0  is  in  the  x-y  plane  referenced  from  the  x-axis.  In  order  to  study  the  simulation  of 
waves  propagating  outward  in  the  z-direction,  the  dispersion  relation  is  rewritten  in  terms 
of  kz,  as  follows: 

K=^-K~k2y  (3.7) 

Simplifying  Equation  (3.7)  gives 

kz=J l-sin2(0)  (3.8) 

It  is  the  approximation  of  the  square  root  in  Equation  (3.8)  that  determines  the  order  of 
the  Mur  operator.  The  radical  is  expanded  in  a  series  using  the  small  argument 
approximation,  thus  limiting  the  angle  of  incidence  6  to  small  values,  and  the  number  of 
terms  kept  in  the  series  yields  the  order  of  the  Mur  operator.  The  first-order  Mur  keeps 
one  term,  leading  to  a  simple  expression  for  kz. 

kz=k  (3.9) 

The  second-order  Mur  keeps  two  terms  and  gives  the  following  approximation  for  kz : 

*,=^l-isin2(0)j  (3.10) 

Recall  that  the  goal  is  to  absorb  a  guided  wave  normally  incident  on  a  z-plane  boundary. 
A  guided  wave  traveling  in  a  dielectric  medium  will  travel  with  a  velocity  slower  than  the 
speed  of  light  in  free  space,  vz  <  ca .  Since  vz  and  kz  are  inversely  related,  a  decrease  in 
vz  means  an  increase  in  kz.  Thus,  kz  should  be  greater  than  k.  This  is  impossible,  given 
Equation  (3.10),  and  explains  the  inability  of  the  second-order  Mur  to  absorb  guided 
waves.  This  serious  limitation  is  not  a  problem  for  the  dispersive  boundary  condition. 
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3.3  Dispersive  Boundary  Condition 

The  dispersive  boundary  condition  analyzed  in  this  work  is  of  the  same  form  as 
the  absorbing  boundary  condition  developed  by  Higdon,  [3.3]  and  [3.4].  The  Higdon 
absorbing  boundary  condition  was  developed  for  radiating  problems,  and  the  parameters 
of  the  ABC  that  can  be  varied  are  the  speeds  of  the  incident  waves  divided  by  the  angles 
of  incidence,  c/cos(0).  This  same  boundary  condition  can  be  used  for  guided  wave 
problems  if  one  interprets  the  wave  speeds  as  normal  velocities  and  modifies  them 
accordingly  [3.5].  The  dispersive  boundary  condition  (DBC)  is  written  as  follows: 


d_  l_d_ 
dz  v,  dt 


A 


L+ll 

dz  v2  dt 


EL  =° 


(3.11) 


The  DBC  has  the  form  of  the  product  of  two  first-order  Mur  operators.  The  differencing 
of  this  expression  is  straightforward  and  results  in  the  following  update  equation  on  a 
uniform  grid: 

pn  _  2F"_i  —  Fn~2  +(v  +  v  V Fn~]  —  Fn  —  Fn~ 2  +  Fn_1  ) 

CJf  ^tf-I  ^M-2  T  V/  1  T  ^M-\  ^M-\  T  ^M-2  ) 

-r,n(E»J-2E«-',+£L)  (3-12) 

where  the  superscripts  represent  the  indexing  in  time,  the  subscripts  represent  the 
indexing  in  space,  and 

\7  —  V  At 

(3.13) 


A z  -  v,A  t 
Yi~  A z  +  V'At 


When  the  operator  (3.11)  is  discretized  on  a  nonuniform  grid,  the  order  in  which  the 
operators  are  applied  has  an  effect  on  the  resulting  update  equation  for  the  boundary. 
Applying  the  operator  with  the  v,  term  first  yields  the  update  equation 

Eu  =2EM_y  —  Em_2  +  (y„  +/l2)(^W  ~  Em-\)  +  {Yu  +Y22){Em-2  ~  ^M-l) 

— YwYni^M  ~  Em-\)~  YuY22(^M-2  —  (3-14) 

A z,  -  vAt 
YiJ~A Zl  +  Vjto 


(3.15) 


where  Az,  is  the  discretization  for  the  cell  closest  to  the  boundary.  Applying  the  operator 
with  v2  first  results  in 

K = -  E«-\ + (r.i + y,2X£« 1  -  ) + (rB + ra,  X^  -  ) 

~YnYn(Eu  ~  Eu-i)~  YnYn[EM-2  ~  ^m-i)  (3-16) 

These  equations  differ  only  in  the  second  and  fourth  multipliers.  Since  it  is  not  obvious 
which  equation  should  perform  better,  a  third  equation  was  written  which  averages  the 
second  and  fourth  multipliers  from  Equations  (3.14)  and  (3.16). 

EM  =2EM_]—EM_2  +  (yu+yl2)(EM  ~  EM_^  —  yuyn{Eu  —  EM_^ 

+[y"  +  ^  ^  ^  +  ^  )(g«-2  _  g»-l) 

+  (3.i7) 

Since  the  multipliers  are  precomputed,  each  of  these  update  equations  requires  only  nine 
additions  and  five  multiplications. 

Numerical  examples  testing  the  performance  of  these  equations  follow  in  Sections 
3.7  and  3.8. 


3.4  The  Modified  Dispersive  Boundary  Condition 

Betz  and  Mittra  [3.6]  introduced  a  damping  factor  into  one  of  the  operators  of  the 
DBC  in  order  to  improve  the  dc  offset  problem  encountered  in  his  research. 
Additionally,  the  damping  factor  can  be  used  in  the  absence  of  the  time  derivative  to 
absorb  evanescent  waves.  The  damping  factor  also  stabilizes  the  DBC  as  is  reported  in 
Section  3.5.  In  this  section,  the  chief  concern  is  the  discretization  of  the  modified 
dispersive  boundary  condition  on  a  nonuniform  grid. 


d  1  d 
—  +  — — +  a 
az  vxdt 


— +  — — )e| 

dz  v2dt)  lz=z' 


=  0 


(3.18) 


Since  evanescent  waves  are  not  considered  in  this  work,  only  one  damping  term  is  used. 
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Discretizing  Equation  (3.18),  assuming  that  the  operator  with  the  damping  factor 
is  applied  first,  results  in  the  following  update  equation 

K  =  £»-,  +  P(K\  -  ) + (r„  +  rniK' 1  -  ■ £«-. ) 

"K/ll  "*■  PY2i){Em-2  ^M-\  )  Y\  \Y  12  (^A/  "**  ) 

-r,  Me"“-2  (3-19) 

where  the  superscripts  represent  time,  the  subscripts  represent  space,  and  y().  is  defined  as 
in  Equation  (3.15).  Following  the  approach  taken  in  [3.7],  the  term  {}  is  defined  as 


P  = 


Az,  +  v,A? 

Az,  +  v,  At(l  +  aAz, ) 


(3.20) 


This  definition  is  obtained  by  applying  the  damping  term  a  only  to  the  term  EnM.  An 
alternative  definition  for  can  be  found  by  applying  central  differencing  in  space  and 
time  to  each  component  of  (3.18)  when  discretizing.  Following  this  procedure  leads  to 
the  definition  in  Equation  (3.21). 


„  AZ]  +  v1Ar(l  -  oAz, ) 
Az,  +  v,  Ar(l  +  aAz, ) 


(3.21) 


There  are  three  parameters  that  can  be  varied  in  Equation  (3.18).  The  effects  of  these 
parameters,  as  well  as  the  definitions  in  Equations  (3.20)  and  (3.21),  on  the  performance 
of  the  DBC  will  be  presented  in  Sections  3.5  -  3.8. 


3.5  Stability 

It  has  been  reported  that  stability  problems  sometimes  arise  when  using 
second-order  and  higher -order  dispersive  boundary  conditions  [3.7],  i.e.,  the  ABC  is  not 
always  reliable.  The  results  of  this  work  also  confirm  this.  Various  explanations  are 
given  for  the  stability  problem.  One  explanation  states  that  the  generation  of  spurious  dc 
signals  depends  on  subtle  factors  such  as  the  exact  order  of  the  additions  in  the  update 
equations,  the  spatial  distribution  of  the  incident  field,  and  the  other  boundary  conditions 
truncating  the  mesh  [3.8].  A  rigorous  approach  to  the  stability  issue  has  been  reported 


regarding  the  stability  of  Liao's  ABC  [3.9].  In  [3.9]  it  has  been  shown  that  one  of  the 
poles  of  the  reflection  coefficient  due  to  the  ABC  is  on  the  unit  circle  in  the  complex 
k-plane,  rendering  the  ABC  marginally  unstable.  It  is  demonstrated  that  perturbing  one 
of  the  multipliers  in  the  differencing  scheme  for  the  Liao  ABC  has  drastic  effects  on  the 
stability  of  the  ABC.  Since  the  DBC  and  Liao's  ABC  are  both  of  the  same  form  as  the 
Higdon  ABC,  this  analysis  can  also  be  applied  to  the  DBC.  Moreover,  it  is  stated  that 
perturbing  the  multiplier  in  order  to  stabilize  the  Liao  ABC  is  equivalent  to  adding  a 
damping  factor  to  the  DBC.  This  will  be  demonstrated  via  a  numerical  example. 

The  mesh  used  to  test  the  stability  of  the  DBC  was  40  by  40  by  41  cells,  and  was 
discretized  with  Ax  =  Ay  =  A z  —  3.750  mm.  The  center  of  the  mesh  was  excited  with  a 
z-directed  point  source  with  the  following  time  variation: 


y(t)  = 


T exp(2)  d 
16  It 


(3.22) 


where 


0.152 


(3.23) 


where  fc  was  chosen  to  be  1.75  GHz  and  At  was  chosen  to  be  6.25  ps.  The  z -component 
of  the  electric  field  was  monitored  15  cells  from  the  source  in  the  y-direction  for  time 
step  0  to  time  step  10000. 


Four  different  ABCs  were  used  in  this  stability  test,  a  first-order  Mur,  and  three 
DBCs  with  8  equal  to  0.00, 0.05,  and  0.10,  where  8  is  defined  in  Equation  (3.24). 

8  =  ocAz  (3.24) 

The  z-component  of  the  electric  field  for  time  steps  0  to  500  is  shown  in  Fig.  3.1. 
Clearly,  the  second-order  DBC  does  a  much  better  job  of  absorbing  the  incident  pulse 
than  the  first-order  Mur.  Figure  3.2  shows  the  electric  field  as  a  function  of  time  for  time 
steps  1500  to  2500.  Note  that  the  electric  field  due  to  the  DBC  without  the  damping  term 
is  growing,  while  the  Mur  and  the  damped  DBC  remain  stable. 


Electric  Field  (V/m) 


Figure  3.1 


Figure  3.2 
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Electric  field  response  near  boundary  for  first-order  Mur  ABC  and  DBC 
with  a  =  8,  for  8  =  0,  8  =  0.05,  and  8  =  0.10. 


Time  Step 


Electric  field  response  near  boundary  for  first-order  Mur  ABC  and  DBC 
with  a  =  8,  for  5=0  (unstable)  and  8  =  0.05  (stable). 


Figure  3.3  demonstrates  the  stability  of  the  damped  DBC,  showing  the  electric 
field  from  time  step  0  to  time  step  10000.  Finally,  Fig.  3.4  shows  that  there  is  a  slight  dc 
offset  when  the  value  of  8  is  not  large  enough.  However,  when  8=  0.10,  there  is  no 
undesirable  dc  offset  and  the  ABC  is  stable. 

The  same  test  was  performed  on  a  nonuniform  mesh,  and  the  DBC  was  found  to 
be  stable  without  the  addition  of  a  damping  term.  Figure  3.5  shows  the  electric  field  5 
cells  from  the  boundary  as  a  function  of  time  for  the  DBC  with  8  equal  to  zero  and  8 
equal  to  0.10.  Although  the  DBC  does  not  become  unstable,  there  is  a  slight  dc  offset 
without  the  damping  term,  as  shown  in  Fig.  3.6.  The  addition  of  the  damping  term 
eliminates  the  dc  offset  as  in  the  uniform  case. 


Time  Step 


Figure  3.3  Electric  field  response  near  the  boundary  for  the  DBC  with  a  =  8,  for 
8  =  0.05  and  8  =  0.10.  The  damping  term  stabilizes  the  DBC. 
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Figure  3.4  Electric  field  response  near  the  boundary  for  the  DBC  with  a  =  <5,  for 
<5  =  0.05  and  8  =  0.10.  The  larger  damping  term  eliminates  the  dc  offset. 
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Figure  3.5  Electric  field  response  near  the  boundary  for  the  DBC  on  a  nonuniform  grid 
with  a  =  8,  for  8  =  0.00  and  <5  =  0.10.  Both  cases  are  stable. 
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Figure  3.6  Electric  field  response  near  the  boundary  for  the  DBC  on  a  nonuniform  grid 
with  a  =8,  for  o  =  0.00  and  8  =  0.10.  The  damping  term  eliminates  the  dc 
offset. 


3.6  Microstrip  Line  with  Uniform  Mesh 

In  this  section,  the  dispersive  boundary  condition  (DBC)  is  tested  on  a  microstrip 
line.  The  width  of  the  line  is  equal  to  the  thickness  of  the  dielectric  substrate,  which  is 
2.54  mm,  as  shown  in  Fig.  3.7.  The  dielectric  constant  of  the  substrate  is  10.2  The  DBC 
is  used  to  truncate  the  end  walls  of  the  microstrip  line,  the  bottom  wall  is  a  perfect  electric 
conductor,  and  the  side  and  top  walls  are  truncated  with  the  first-order  Mur  ABC.  The 
mesh  is  discretized  with  Ax  =  Ay  =  Az  =  0.3715  mm,  and  the  time  step  At  is  0.6053  ps, 
which  is  0.99  times  the  limiting  value  from  the  Courant  stability  criterion.  The  structure 
is  excited  with  a  window  function  with  a  3  dB  cutoff  frequency  of  8  GHz.  To  test  the 
effectiveness  of  the  DBC,  the  problem  is  run  with  two  meshes.  The  test  mesh  is  24  by  44 
by  270  and  the  larger  reference  mesh  is  24  by  44  by  500.  The  voltage  is  monitored  along 
the  line  as  a  function  of  time.  The  solution  on  the  reference  line  is  considered  to  be  exact 
since  the  simulation  is  stopped  before  reflections  from  the  far  end  of  the  line  can  reach 
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the  point  where  the  voltage  is  monitored  and  corrupt  the  signal.  The  difference  between 
the  reference  solution  and  the  test  solution  will  show  the  effectiveness  of  the  ABC  being 
tested. 

The  normalized  reflected  voltage  is  plotted  as  a  function  of  time  in  Fig.  3.8.  The 
curve  labeled  "mur"  is  the  result  using  a  first-order  Mur  ABC  on  the  test  boundary.  The 
ABC  is  implemented  so  that  the  speed  of  the  wave  is  determined  by  the  material  property 
of  the  medium  in  a  given  cell.  For  example,  cells  in  air  are  updated  with  the  velocity  of 
free  space,  whereas  cells  in  the  dielectric  are  updated  with  a  velocity  of  c/ -JfTr.  The  peak 
of  the  normalized  reflection  in  the  time  domain  is  around  0.025.  The  curve  labeled 
"dbcO"  is  the  result  of  setting  <5  =  0,  eri  =  7.12,  and  eh  =  8.5.  For  "dbcl",  8  =  0.1, 
erj  =  7.12,  and  eh  -  8.5,  and  for  "dbc2",  8  =  0.1,  en  =  8.5,  and  e h  =  7.12.  All  three  of 
the  DBC  results  appear  to  be  considerable  improvements  when  compared  with  the  Mur 
result. 


Figure  3.7  Geometry  of  the  microstrip  line  used  to  test  the  DBC.  The  relevant 
dimensions  are  labeled  in  the  figure. 
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Time  Step 

Figure  3.8  Normalized  reflected  voltage  on  the  microstrip  line  as  a  function  of  time: 

"mur"  is  first-order  Mur;  "dbcO"  has  5=0,  erl  =  7.12,  and  er2  =  8.5;  "dbcl" 
has  8  =  0.1,  £r]  =  7.12,  and  er2  =  8.5;  "dbc2"  has  8  =  0.1,  en  =8.5,  and 
Cr2  =7.12. 

Another  useful  way  to  ascertain  the  effectiveness  of  the  ABCs  under  test  is  to 
define  a  reflection  coefficient  due  to  the  ABC.  The  test  voltage  is  subtracted  from  the 
reference  voltage,  giving  the  reflected  voltage.  The  test  voltage  and  reflected  voltages  are 
transformed  to  the  frequency  domain,  and  the  reflection  coefficient  is  then  defined  as  the 
ratio  of  the  reflected  voltage  to  the  incident  voltage.  The  magnitude  of  the  reflection 
coefficient  is  plotted  as  a  function  of  frequency  in  Fig.  3.9,  where  it  can  be  seen  that  the 
best  performing  ABC  is  the  DBC  with  8  set  to  zero.  Although  this  yields  the  best  result, 
it  can  potentially  lead  to  late  time  instabilities,  as  was  demonstrated  in  Section  3.5.  Since 
the  8  term  is  applied  in  the  first  operator,  the  choice  of  ert  and  eh  is  significant.  Clearly, 
it  is  best  to  choose  the  smaller  er  value  for  eri .  This  results  in  a  reflection  coefficient 
with  magnitude  less  than  0.01  over  the  1  GHz  to  16  GHz  range. 
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Figure  3.9  Magnitude  of  reflection  coefficient  as  a  function  of  frequency  for  the 
microstrip  line.  Curves  are  labeled  as  in  Fig.  3.8. 


The  error  due  to  the  ABCs  can  also  be  quantified  using  error  norms.  Two 
different  norms  are  used  in  this  analysis,  an  I,  norm  and  an  norm.  The  norms  are 
defined  as  follows: 


(  n 
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Ki=i 

J 

(3.25) 


where  vr  and  v,  are  the  reference  and  test  voltages,  respectively.  Error  norms  are  listed 
in  Table  3.1.  The  numerical  data  show  the  same  trend  as  the  graphical  data  presented. 


Table  3.1  Error  for  DBC  tested  on  microstrip  line 


DBC 

s 

er2 

LI  (*10~2) 

L2  (*10'4) 

0 

0.0 

7.12 

8.50  ! 

0.392367 

1.765111 

1 

0.1 

7.12 

8.50 

0.457797 

2.068264  4 

2 

0.1 

8.50 

7.12 

0.778930 

3.206607 

First-order  Mur 

1.429872 

6.966470 
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3.7  Microstrip  Line  with  Nonuniform  Mesh 

The  effectiveness  of  the  dispersive  boundary  condition  (DBC)  when  applied  to  a 
uniform  mesh  was  demonstrated  in  previous  sections.  In  this  section,  the  DBC  is  applied 
to  a  microstrip  line  with  a  nonuniform  mesh.  Again,  two  meshes  are  used  to  test  the 
effectiveness  of  the  DBC.  The  large  reference  mesh  is  24  by  44  by  500  cells  and  the  test 
mesh  is  24  by  44  by  270  cells.  The  geometry  is  excited  and  the  voltage  nonitored  in  the 
same  manner  as  in  the  uniform  case.  In  order  to  isolate  the  effects  of  the  nonuniform 
mesh  and  nonuniform  boundary  condition,  only  the  last  four  cells  in  the  z-direction  are 
nonuniform  and  the  growth  rate  between  neighboring  cells  is  at  most  1.20.  Three 

parameters  are  varied  when  testing  the  performance  of  the  DBC.  These  parameters  are 
ef) ,  er2 ,  and  S,  where  S  is  defined  in  Equation  (3.24)  and  a  appears  with  en  as  shown 

in  Equation  (3.18). 

Different  values  and  combinations  of  er  and  S  were  tested,  and  time  domain 
results  are  plotted  in  Figs.  3.10  -  3.12.  In  these  figures,  the  difference  between  the 
reference  voltage  and  test  voltage  is  normalized  by  the  maximum  value  of  the  reference 
voltage.  The  values  of  eh  and  eri  are  7.12  and  8.5,  respectively,  for  cases  0,  1,  3,  and  4. 
For  case  2,  e  =  8.5  and  eh  -  7.12.  The  <5  values  for  cases  0  through  4  are  0.0,  0.10, 
0.10,  0.05,  and  0.01,  respectively.  Figure  3.10  shows  that  case  0  produces  less  reflection 
than  case  1.  The  result  from  the  first-order  Mur  ABC  is  also  provided  in  Fig.  3.10.  The 
peak  in  the  time  domain  reflection  from  the  first-order  Mur  ABC  is  more  than  three  times 
greater  than  the  largest  peak  from  the  DBC.  Figure  3.1 1  compares  cases  2  through  4  and 
clearly  shows  that  case  2  yields  poor  results.  Comparing  the  results  of  Fig.  3.10  with  the 
results  in  Fig.  3.11  shows  that  although  the  reflection  from  case  2  is  large,  it  is  still  an 
improvement  over  the  first-order  Mur  result.  The  time  range  is  narrowed  in  Fig.  3.12  and 
the  effect  of  varying  S  is  shown.  The  smaller  values  of  S  produce  smaller  reflections  in 
this  time  range. 
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Figure  3.10  Normalized  reflected  voltage  as  a  function  of  time  for  the  microstrip  line  on 
a  nonuniform  grid.  The  first-order  Mur  result  is  labeled  "mur,"  e  =7.12, 
er2  =  8.5,  "dbcO"  has  8  =  0,  and  "dbcl"  has  S  =  0.1 
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Figure  3.11  Normalized  reflected  voltage  as  a  function  of  time  for  the  microstrip  line  on 
a  nonuniform  grid.  "dbc2"  has  S  =  0.1,  erl  =  8.5,  and  er2  =  7.12,  "dbc3" 
has  S  =  0.05,  erl  =  7.12,  and  er2  =  8.5,  "dbc4"  has  S  =  0.01,  eH  =  7.12,  and 
er2  =  8.5. 
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Figure  3.12  Normalized  reflected  voltage  as  a  function  of  time  for  the  microstrip  line  on 
a  nonuniform  grid.  Curves  labeled  as  in  Figs.  3.10  and  3.11. 

Frequency  domain  plots  of  the  magnitude  of  the  reflection  coefficient  are  shown 
in  Figs.  3.13  and  3.14.  Figure  3.13  shows  that  when  er>  =  7.12  and  efj  =  8.5,  the  DBC 
outperforms  the  first-order  Mur  by  a  factor  of  around  5  or  6  in  the  frequency  range  of 
4  GHz  to  10  GHz.  When  er_  =  8.5,  eri  =  7.12,  and  8  =  0.10,  the  DBC  is  still  much  more 
effective  than  the  first-order  Mur,  which  is  not  obvious  from  examining  the  time  domain 
results.  Figure  3.14  compares  the  effects  of  varying  S.  Again,  the  smaller  8  produces 
smaller  reflections. 

As  with  the  uniform  mesh  microstrip  line  test,  the  errors  can  be  quantified  with 
error  norms.  The  numerical  data  presented  in  Table  3.2  show  the  same  trends  as  the 
graphical  data.  To  minimize  the  reflection  from  the  absorbing  boundary,  the  smaller 
value  of  the  two  dielectric  constants  should  be  used  in  the  first  operator,  and  the  damping 
term  should  be  small,  but  remain  nonzero  in  order  to  prevent  late-time  instabilites. 


Figure  3.13  Magnitude  of  reflection  coefficient  as  a  function  of  frequency  for  the 
microstrip  line  on  a  nonuniform  grid.  Curves  labeled  as  in  Figs.  3.10  and 
3.11. 


Figure  3.14  Magnitude  of  reflection  coefficient  as  a  function  of  frequency  for  the 
microstrip  line  on  a  nonuniform  grid.  Curves  labeled  as  in  Figs.  3.10  and 
3.11. 


Table  3.2  Error  for  DBC  tested  on  microstrip  line  with  nonuniform  mesh 


DBC 

6 

^rl 

er2 

LI  (*10~2) 

L2  (*10*4) 

0 

0.00 

7.12 

8.50 

0.452295 

1.448510 

1 

0.10 

7.12 

8.50 

0.491422 

1.997053 

2 

0.10 

8.50 

7.12 

0.814429 

3.536524 

3 

0.05 

7.12 

8.50 

0.431705 

1.701391 

4 

0.01 

7.12 

8.50 

0.453550 

1.488936 

First-order  Mur 

1.473694 

7.053369 

3.8  Angle  Absorbing  Boundary  Condition 

The  form  of  the  dispersive  boundary  condition  (DBC)  allows  one  to  choose  the 
velocity  which  will  optimize  the  performance  of  the  boundary  condition  for  any  given 
problem.  This  is  important  since  it  can  then  be  used  for  both  guided  wave  and  radiating 
wave  problems.  When  the  DBC  is  used  for  a  radiation  problem,  it  will  be  referred  to  as 
the  angle  absorbing  boundary  condition  (AABC). 

In  order  to  test  the  AABC,  a  simple  point  source  radiating  in  free  space  was 
analyzed.  A  large  mesh  with  104  by  104  by  103  cells  provided  the  reference  solution  and 
a  smaller  mesh  with  54  by  54  by  53  cells  was  used  to  test  the  AABC.  The  majority  of  the 
test  mesh  is  uniform  with  Ax  =  Ay  =  Az  =  9.000  mm.  The  outer  four  layers  in  each 
direction  are  nonuniform  with  a  maximum  growth  rate  of  1.2.  The  time  step  is  chosen  to 
be  15.625  ps  and  the  center  of  the  mesh  is  excited  with  the  following  z-directed  source: 

Ez(n)  =  A  *  (10  -  15cos(<p)  +  6cos(2p)  -  3cos(3<p))  (3.26) 

where  A  =  1.0/320.0,  (p  =  27tfrnAt,  and  fr  =  1  GHz.  The  source,  shown  in  Fig.  3.15,  is 
nonzero  for  64  time  steps. 
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Figure  3.15  Time  signature  of  the  point  source  excitation  used  to  test  the  angle  absorbing 
boundary  condition. 


Error  analysis  is  carried  out  as  a  function  of  time  and  position.  The  z-component 
of  the  electric  field  is  monitored  from  time  step  80  to  time  step  160  on  the  plane  two  cells 
from  the  y  =  0  boundary.  This  allows  the  source  pulse  to  pass  through  the  boundary  and 
shows  the  effect  of  the  imperfect  radiating  boundary  condition.  When  the  error  is  shown 
as  a  function  of  time,  the  absolute  value  of  the  difference  between  the  reference  solution 
and  the  test  solution  is  summed  at  each  location  on  the  test  plane  and  normalized  by  the 
maximum  value  of  the  sum  of  the  magnitudes  of  the  field  values  in  the  reference  solution 
in  the  time  range  of  interest.  The  expression  for  the  error  is  given  in  Equation  (3.27) 


Error(t )  = 


maxi 


(3.27) 


where  vr  is  the  reference  field  value,  v,  is  the  test  field  value,  i  and  j  are  spatial  indices, 
and  t  is  the  time  index.  When  the  error  is  presented  as  a  function  of  position,  the  sum  of 
the  absolute  value  of  the  difference  between  the  reference  solution  and  the  test  solution 


from  time  step  80  to  time  step  160  is  normalized  by  the  sum  of  the  absolute  value  of  the 
reference  solution  over  that  same  time  interval  for  a  given  position,  as  shown  in 
Equation  (3.28). 


Error(p) = 


ZJK  M-y,M 


(3.28) 


where  p  is  the  position  index.  The  position  moves  along  a  diagonal  from  the  lower  left 
comer  to  the  center  of  the  test  plane.  This  corresponds  to  moving  from  position  0  to 
position  25  as  shown  in  Fig.  3.16.  A  diagonal  line  was  chosen  in  order  to  gauge  the 
performance  of  the  ABCs  being  tested  as  a  function  of  position.  This  is  important 
because  the  performance  of  ABCs  in  comers  is  of  interest.  Using  monitoring  locations 
along  a  diagonal  gives  a  wider  range  of  incident  angles  upon  the  boundary  as  opposed  to 
a  line  parallel  to  one  of  the  edges  of  the  mesh. 


Figure  3.16  The  geometry  used  to  test  the  Angle  Absorbing  Boundary  Condition  is  a 
point  source  radiating  in  free  space.  The  fields  are  monitored  on  a  line  in 
the  y  =  2  plane. 
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As  discussed  previously,  the  order  of  the  implementation  of  the  operators  does 
make  a  difference  on  the  nonuniform  grid,  whether  or  not  the  damping  term  is  present. 
The  AABC  was  tested  with  no  damping  term  using  two  different  angles,  one  of  which 
was  always  zero.  For  a  given  angle,  there  were  three  possible  implementations.  One 
with  v,  modified  by  the  angle  and  v2  equal  to  c,  one  with  v2  modified  by  the  angle,  and 
one  with  the  multipliers  determined  by  averaging  the  multipliers  from  the  first  two  cases. 
The  nonzero  angles  tested  were  15°,  30°,  and  45°.  The  results  are  compared  with  the 
first-order  Mur. 

Results  for  normalized  error  as  a  function  of  time  are  shown  in  Figs.  3.17  -  3.20. 
Figure  3.17  shows  results  in  which  the  v,  term  is  modified  by  the  angle  as  shown  on  the 
plot  legend  and  "mur"  refers  to  the  first-order  Mur  ABC.  Figures  3.18-  3.20  compare  the 
errors  obtained  using  the  three  different  implementations  of  the  AABC  for  a  given  angle. 
The  best  results  are  achieved  when  v,  is  modified  by  a  45  degree  angle. 


Figure  3.17  Normalized  error  as  a  function  of  time,  "mur"  is  first-order  Mur,  d2  =  0,  and 
0,  varies  in  degrees  as  shown  in  the  legend. 
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Figure  3.18  Normalized  error  as  a  function  of  time,  "vl-15"  has  0,  =  15°  and  02  = 
"v2-15"  has  0,  =  0°  and  02  =  15°;  "a-15"  averages  the  multipliers. 


Figure  3.19  Normalized  error  as  a  function  of  time,  "vl-30"  has  0,  =  30°  and  02  = 
"v2-30"  has  0,  =  0°  and  02  =  30°;  "a-30"  averages  the  multipliers. 
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Figure  3.20  Normalized  error  as  a  function  of  time,  "vl-45"  has  0,  =  45°  and  d2  =  0°; 
"v2-45"  has  6X  =  0°  and  02  =  45°;  "a-45"  averages  the  multipliers. 


Monitoring  the  error  as  a  function  of  position  shows  the  benefits  of  using  the 
AABC.  Figures  3.21  -  3.24  show  the  normalized  error  as  a  function  of  position  for 
various  angles.  The  optimal  results  are  obtained  when  an  angle  of  45  degrees  is  used  in 
conjunction  with  the  v,  operator.  This  is  particularly  apparent  near  the  comer  of  the 
mesh,  where  the  45  degree  angle  causes  the  normalized  error  to  be  less  than  5%.  The 
error  from  the  first-order  Mur,  on  the  other  hand,  is  on  the  order  of  33%  near  the  comer. 
This  is  a  key  point  because  the  second-order  Mur  cannot  be  implemented  directly  in  the 
comers,  and  when  an  averaged  first-order  Mur  is  used,  significant  error  is  intoduced  into 
the  numerical  data.  As  expected,  the  performance  near  the  comer  improves  as  the  angle 
varies  from  zero  to  45  degrees.  It  is  interesting  to  note  that  averaging  the  multipliers 
gives  better  results  than  modifying  the  v2  operator  only  when  the  45  degree  angle  is  used. 
This  can  be  seen  from  Figs.  3.22  -  3.24. 


Normalized  Error  £  Normalized  Error 
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The  total  error  can  be  measured  in  terms  of  the  L,  norm  as  defined  previously. 
By  summing  the  error  as  a  function  of  both  position  and  time,  a  single  numerical  value  is 
obtained.  This  procedure  was  followed  in  order  to  quantify  the  total  error  for  the 
different  implementations  and  parameter  values.  The  calculated  norms  are  normalized  by 
the  I,  norm  due  to  the  first-order  Mur  ABC.  Figure  3.25  shows  the  normalized  total 
error  as  a  function  of  angle.  The  v,  operator  has  a  45°  angle  and  the  angle  in  the  v2 
operator  varies.  Clearly  it  is  best  to  have  0,  =  45°  and  02  =  0°  to  minimize  the  total  error. 


Figure  3.25  Error  normalized  to  error  due  to  first-order  Mur.  0,  =  45°  and  02  varies. 

The  total  error  can  also  be  plotted  as  a  function  of  the  damping  factor,  8.  Two 
possible  differencing  schemes  were  presented  in  a  previous  section.  These 
implementations  of  the  modified  DBC  will  be  referred  to  as  "dl"  and  "d2,"  where  "dl"  is 
given  by  Equation  (3.20)  and  "d2"  is  given  by  Equation  (3.21).  The  L,  norms  of  total 
error  are  normalized  by  the  total  error  due  to  a  first-order  Mur  ABC  and  plotted  as  a 
function  of  8  in  Fig.  3.26.  As  expected,  the  "d2"  scheme  is  more  sensitive  to  the  value  of 
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8  since  it  appears  in  both  the  numerator  and  denominator  of  the  multipliers.  In  the  "dl" 
scheme,  the  damping  factor  is  found  only  in  the  denominator;  hence,  the  error  is  less 
sensitive  to  variations  of  8.  Figure  3.26  shows  that  error  is  minimized  for  "dl"  when  8 
is  in  the  range  of  0.15  to  0.20.  Error  is  minimized  for  "d2"  when  8  is  around  0.08.  The 
performance  of  the  "dl"  implementation  is  only  slightly  better  than  for  the  "d2"  scheme. 

The  addition  of  the  damping  term  not  only  stabilizes  the  AABC,  but  also 
improves  the  performance.  Figures  3.27  and  3.28  show  the  normalized  error  as  a 
function  of  position  for  the  "d2"  implementation,  with  0,  =  45°  and  02  =  0°.  Three 
values  of  8  are  used  as  shown  on  the  plot  legend.  Figure  3.27  displays  the  L,  norm  error 
normalized  by  the  sum  of  the  reference  solution,  while  Fig.  3.28  shows  the  A,  norm  error 
normalized  by  the  square  root  of  the  sum  of  the  squares  of  the  field  values  in  the 
reference  solution.  These  curves  show  the  same  type  of  behavior  as  the  curves  in  Figs. 
3.21  -  3.24.  One  interesting  feature  of  these  results  is  that  the  damping  term  reduces  the 
error  for  positions  6  through  17. 


Figure  3.26  Error  normalized  to  error  due  to  first-order  Mur.  0,  =  45°  and  02  =  0°.  "dl" 
uses  Equation  (3.20)  and  "d2"  uses  Equation  (3.21). 
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Figure  3.27  L,  normalized  error  as  a  function  of  position.  =  45°  and  d2  =  0°,  "mur" 
is  first-order  Mur;  AABC  has  8  values  as  shown  in  the  legend. 
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Figure  3.28  L2  normalized  error  as  a  function  of  position.  0,  =  45°  and  d2  =  0°,  "mur" 
is  first-order  Mur;  AABC  has  <5  values  as  shown  in  the  legend. 
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The  improvement  is  also  apparent  when  looking  at  the  error  as  a  function  of  time. 
In  Fig.  3.29,  the  normalized  error  from  8  equal  to  zero  is  compared  with  results  when  8 
equals  0.10  for  two  cases.  In  case  1,0,=  45°  and  02  =  0°  and  in  case  2,  0,  =  0°  and 
02  =  45°.  Case  1  is  labeled  as  vl-45  on  the  plot  legend,  and  the  "d"  indicates  the  nonzero 
S  value.  The  damping  term  actually  makes  the  performance  worse  for  case  2,  but  case  2 
is  inferior  to  case  1  in  performance.  Case  1  is  improved  with  the  damping  factor,  as 
desired.  Figure  3.30  shows  the  effect  of  varying  the  value  of  S  for  0,  =  45°  and  02  =  0°. 
The  normalized  error  is  plotted  as  a  function  of  time;  it  appears  that  the  best  results  occur 
when  8  is  around  0.20.  This  agrees  with  the  results  in  Fig.  3.26. 

Error  norms  were  calculated  and  are  presented  in  Table  3.3.  These  values  agree 
with  the  graphical  data,  showing  that  the  best  performance  is  obtained  when  using 
0,  =  45°,  02  =  0°,  and  a  nonzero  8  which  can  be  chosen  according  to  the  results  in 
Fig.  3.26. 
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Figure  3.29  Normalized  error  as  a  function  of  time,  "vl-45"  indicates  that  0,  =  45°  and 
02  =  0°  and  "d"  means  that  8  =  0.10. 
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Figure  3.30  Normalized  error  as  a  function  of  time.  0,  =  45°,  02  =  0°  and  8  varies  as 
shown  in  the  legend. 


Table  3.3  Error  for  AABC  for  a  point  source  radiating  in  free  space.  Error  is  normalized 

to  the  error  from  the  first-order  Mur  ABC. 


8 

0,  (degrees) 

02  (degrees) 

Error 

0.1 

45 

0 

0.154983 

0.0 

45 

0 

0.205506 

0.1 

30 

0 

0.223243 

0.0 

45,0 

average 

0.253852 

0.0 

30 

0 

0.258989 

0.1 

15 

0 

0.268537 

0.1 

0 

30 

0.276025 

0.1 

0 

15 

6.277707 

0.0 

0 

30 

0.282492 

0.0 

15 

0 

0.298200 

0.0 

0 

15 

0.301365 

0.0 

0 

45 

0.311096 

“0.0 

15,0 

average 

0.324217 

0.1 

0 

45 

0.344036 

0.0 

30,0 

average 

0.368374 

First-order  Mur 

1.000000 
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3.9  Conclusions 

The  study  of  the  Dispersive  Boundary  Condition  was  motivated  by  a  discussion  of 
the  shortcomings  of  the  Mur  Absorbing  Boundary  Conditions.  The  stability  issue  was 
discussed  in  conjunction  with  the  Modified  Dispersive  Boundary  Condition.  It  was 
demonstrated  that  the  DBC  can  be  used  to  effectively  absorb  guided  waves  and  radiating 
waves.  The  DBC  was  applied  to  nonuniform  meshes  and  tested  with  a  microstrip  line 
problem  and  with  a  radiating  point  source  problem.  Several  different  implementations  of 
the  DBC  on  a  nonuniform  mesh  were  presented.  The  parameters  of  the  DBC  were  varied 
and  results  were  presented  showing  the  optimal  perfomance  of  the  DBC  in  terms  of 
minimizing  reflection  error.  Several  different  forms  of  error  analysis  were  presented 
enabling  one  to  gauge  the  performance  of  the  DBC  in  terms  of  the  error  as  a  function  of 
both  time  and  position. 
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CHAPTER  4 


CHARACTERIZATION  OF  MICROSTRIP  LINES  AND  DISCONTINUITIES 
USING  NONUNIFORM  ORTHOGONAL  FDTD 


4.1  Introduction 

Accurate  characterization  of  frequency -dependent  transmission  lines  and 
transmission  systems  is  an  important  topic  in  the  area  of  millimeter-wave,  microwave, 
and  digital  circuit  design,  in  which  it  is  desirable  to  minimize  the  size  of  circuits  and 
electronic  packages.  This  concern  with  size  reduction  often  forces  lines  to  bend  and  to  be 
close  to  one  another,  which  introduces  frequency-dependent  effects.  Moreover,  as  clock 
rise  times  become  smaller,  the  electromagnetic  effects  of  the  wave  propagating  along  the 
signal  lines  and  throughout  the  circuit  become  increasingly  important.  The  FDTD 
method  is  an  excellent  tool  for  solving  these  types  of  problems  because  it  provides  a 
full -wave  solution  for  a  wide  range  of  frequencies  with  only  a  single  simulation.  Circuit 
simulators,  although  very  fast,  on  the  other  hand,  use  approximations,  which  are  only 
valid  under  somewhat  stringent  restrictions  on  frequency  range  and  geometrical  features 
depending  upon  the  circuit  elements.  Many  of  the  approximations  are  quasi-static,  and 
thus  cannot  be  used  for  high  frequencies.  Additionally,  in  many  instances,  circuits  are 
composed  of  complex  transitions  which  cannot  be  modeled  using  circuit  simulators.  This 
is  not  a  limitation  for  the  FDTD  method,  which  can  solve  arbitrarily  complicated 
geometries  including  active  and  passive  loads  [4.1]  limited  only  by  the  memory 
limitations  of  the  computer  used  for  the  FDTD  simulation. 

In  this  chapter,  the  frequency -dependent  characteristics  of  a  uniform  microstrip 
line  are  calculated,  thus  demonstrating  the  applicability  of  the  FDTD  method.  Use  of  a 
Prony  technique  with  the  FDTD  method  is  presented  and  applied  to  two-port  scattering 
parameter  calculations.  A  special  update  equation  for  triangular  metallization  is 


introduced  and  shown  to  increase  the  accuracy  of  the  FDTD  method.  Finally,  the  FDTD 
method  and  Prony  technique  are  combined  to  analyze  a  grounded  coplanar  waveguide  to 
microstrip  line  transition. 

4.2  Uniform  Microstrip  Line 

In  this  section,  the  uniform  microstrip  line  of  Fig.  4. 1  is  analyzed  using  the  FDTD 
method.  The  width  of  line,  W,  is  150  pm,  the  height  of  the  substrate,  H,  is  100  pm,  and 
the  dielectric  constant,  er,  is  13.0.  These  same  line  parameters  are  used  in  [4.2]  and 
[4.3].  A  uniform  cell  discretization  of  Ax  =  Ay  =  Az  =  12.5  pm  was  used  and  At  was 
set  to  20.8333  fs.  The  line  was  excited  with  a  Blackman-Harris  window  function  with  a 
3  dB  cutoff  frequency  of  around  240  GHz,  and  the  voltage  and  current  were  monitored  at 
several  locations  along  the  line.  The  mesh  used  to  describe  the  microstrip  line  was 
30  x  108  x  160  cells.  The  bottom  wall  was  a  plane  of  perfect  electric  conductor,  and  the 
top  and  side  walls  were  truncated  with  the  first-order  Mur  absorbing  boundary  condition 
(ABC).  The  absorbing  boundary  conditions  used  on  the  walls  that  terminate  the  line  were 
varied.  Figure  4.2  shows  the  magnitude  of  the  reflection  coefficient  due  to  the  absorbing 
boundary  as  a  function  of  frequency  for  the  first-order  Mur  ABC  and  two 
implementations  of  the  dispersive  boundary  condition  (DBC). 


Figure  4.1  Geometry  of  the  uniform  microstrip  line.  The  important  parameters  are  the 
line  width,  W,  the  substrate  height,  H,  and  the  dielectric  constant,  er. 
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Figure  4.2  Magnitude  of  reflection  coefficient  of  the  imperfect  absorbing  boundary 
conditions  used  to  truncate  the  mesh  describing  the  line  in  Fig.  4.1.  Mur 
refers  to  a  first-order  Mur.  The  dispersive  boundary  conditions  have  e  =9.4 

and  £  =1 1.2.  DBC1  has  5=0.10  and  DBC2  has  5=0.01. 

r2 

Clearly,  it  is  best  to  use  a  small  value  of  5  when  using  the  DBC  to  minimize  the 
amount  of  reflection.  The  curve  labeled  "DBC2"  in  Fig.  4.2  shows  reflections  of  less 
than  0.5%  over  the  20  GHz  to  160  GHz  frequency  range.  When  compared  with  the 
roughly  3%  error  of  the  first-order  Mur  over  the  same  frequency  range,  the  DBC  provides 
a  significant  reduction  in  the  fields  reflected  from  the  ABC.  Although  the  magnitude  of 
the  error  shown  in  Fig.  4.2  for  the  first-order  Mur  is  small,  it  corrupts  the  solution  enough 
to  cause  oscillations  in  the  effective  dielectric  constant  and  impedance  calculations  [4.4]. 
These  oscillations  can  be  removed  by  using  a  technique  based  on  Prony’s  method  [4.5]. 
This  method  is  presented  in  the  next  section,  and  was  used  to  calculate  the  reflection 
coefficient  shown  in  Fig.  4.2,  as  well  as  the  effective  dielectric  constant  shown  in  Fig.  4.3 
and  the  characteristic  impedance  shown  in  Fig.  4.4.  There  is  no  noticeable  difference 
between  the  values  calculated  using  the  different  ABCs,  thus  demonstrating  the 
robustness  of  the  technique.  The  agreement  with  the  results  of  Becker  is  quite  good. 


Effective  Dielectric  Constant,  e 


Frequency  (GHz) 

Figure  4.3  Effective  dielectric  constant  for  the  microstrip  line  of  Fig.  4.1.  Absorbing 
boundary  parameters  are  as  noted  in  Fig.  4.2.  The  published  results  of 
Becker  [4.3]  are  provided  as  a  reference. 


Figure  4.4  Characteristic  impedance  of  the  microstrip  line  of  Fig.  4.1.  Absorbing 
boundary  parameters  are  as  noted  in  Fig.  4.2.  The  published  results  of 
Becker  [4.3]  are  provided  as  a  reference. 
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4.3  Prony's  Method  for  Scattering  Parameter  Extraction 

The  usual  method  [4.6]  for  determining  scattering  parameters  of  microstrip 
discontinuities  requires  two  FDTD  simulations.  The  first  run  simulates  the  uniform 
microstrip  line  in  order  to  determine  the  incident  voltage  and/or  current.  The  second  run 
simulates  the  line  with  the  discontinuity  and  finds  the  total  voltage  and/or  current  on  the 
source  side  and  the  transmitted  voltage  on  the  other  side.  The  reflected  voltage  is  simply 
the  difference  between  the  total  voltage  and  the  incident  voltage.  The  scattering 
parameters  are  defined  in  terms  of  the  Fourier  transforms  of  the  time  signatures  of  the 
incident,  reflected,  and  transmitted  voltages  as  shown  in  (4.1)  and  (4.2). 
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SAf)  =  YJW  (4-2) 

It  is  possible  to  calculate  the  scattering  parameters  with  only  one  FDTD 
simulation  and  to  eliminate  the  adverse  effects  from  the  imperfect  ABCs  [4.5].  This  is 
achieved  by  using  a  method  based  on  Prony’s  method  and  transmission  line  equations. 
Consider  the  general  miscrostrip  discontinuity  problem  shown  in  Fig.  4.5,  where  the 
voltage  on  the  source  side  of  the  discontinuity  consists  of  a  forward  traveling  wave  and  a 
backward  traveling  wave,  as  does  the  voltage  on  the  other  side  of  the  discontinuity.  If  the 
line  runs  in  the  z-direction  and  the  forward  traveling  source  is  defined  as  A ,  then  the 
voltage  on  Side  Two  can  be  written  as 

V2(z)  =  S214e-*+r,»A.4e-!’tV«  (4.3) 

where  roic  is  the  reflection  due  to  the  imperfection  of  the  absorbing  boundary,  is  the 
length  of  the  line  from  reference  plane  2  to  the  boundary,  y  is  the  propagation  constant, 
and  S21  is  one  of  the  unknown  scattering  parameters.  Similarly,  the  voltage  on  Side  One 
can  be  written  as 
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v,(z)  =  (A+  r .AAe-2^  y*  +  S„Ae **  (4.4) 

where  L,  is  the  distance  from  reference  plane  1  to  the  boundary,  and  Su  is  the  other 
unknown  scattering  parameter.  Since  the  same  numerical  boundary  condition  is  applied 
to  all  boundaries,  and  the  distances  from  the  reference  planes  to  the  boundaries  are  the 
same,  ra6c  is  the  same  on  side  one  as  on  side  two.  It  should  be  noted  that  higher -order 
products  of  scattering  parameters  and  reflections  due  to  the  boundaries  are  being 
neglected. 

Both  Equations  (4.3)  and  (4.4)  are  of  the  form  V(z )  =  A+e~yz  +A~e+r\  where  A+ 
and  A~  represent  the  complex  amplitude  of  the  forward  and  backward  traveling  waves, 
respectively.  These  amplitudes,  together  with  the  propagation  constant,  give  three 
unknowns.  Thus,  if  the  voltage  is  monitored  at  three  locations  along  the  line,  the  three 
unknowns  can  be  determined  in  the  following  manner.  First,  consider  the  voltage  on  Side 
Two.  The  voltage  is  monitored  at  three  z-locations  in  (4.3)  resulting  in  a  system  of 
equations  which  is  solved  using  the  Prony  method  to  determine  y ,  S2lA,  and  Yabc.  Next, 
three  voltages  on  Side  One  with  Tabc  from  Side  Two  are  used  to  determine  A  and  Su. 
Finally,  the  value  of  A  determined  on  Side  One  is  used  to  calculate  S21.  This  method  is 
twice  as  efficient  as  the  usual  method  for  calculating  scattering  parameters.  Numerical 
results  are  presented  in  the  next  section. 
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Figure  4.5  General  microstrip  line  discontinuity  with  scattering  parameter  reference 
planes. 


4.4  Microstrip  Line  with  Right-Angle  Bend 

A  common  microstrip  line  discontinuity  is  the  right-angle  bend.  The  scattering 
parameters  for  this  discontinuity  will  be  calculated  in  this  section.  The  geometry  of  the 
right-angle  bend  is  shown  in  Fig.  4.6.  The  line  is  4.4  mm  wide  and  55  mm  in  length  on 
both  sides  of  the  bend.  The  substrate  has  thickness  1.6  mm  and  er  =  2.62.  These 
dimensions  result  in  an  impedance  of  approximately  50  Q  near  2.5  GHz.  The 
discretization  of  the  mesh  is  Ax  =  0.533  mm,  Ay  =  Az  =  0.55  mm,  and  the  time  step  is 
At  =  907.1201  fs.  The  reference  planes  for  calculating  the  scattering  parameters  are 
located  32  cells  (17.6  mm)  away  from  the  bend.  The  computational  domain  is 
140  x  140  x  15  cells,  and  the  first-order  Mur  ABC  is  used  to  truncate  the  mesh  except  for 
the  x  =  0  boundary  which  is  a  perfect  electric  conductor.  The  microstrip  line  is  excited 
1 1  mm  from  the  z  =  0  boundary  with  a  Blackman-Harris  window  function  with  a  3  dB 
cutoff  frequency  of  10  GHz. 


Figure  4.6  Geometry  of  microstrip  line  with  right-angle  bend  discontinuity. 


The  frequency  domain  results  in  Figs.  4.7,  4.8,  and  4.9  compare  the  scattering 
parameters  obtained  when  analyzing  the  right-angle  bend  using  the  FDTD  method  with 
two  runs,  using  one  FDTD  run  and  the  Prony  technique,  and  using  the  circuit  simulator 
Touchstone  by  EEsof  [4.7].  Results  from  Touchstone  are  limited  to  14  GHz,  and  they  do 
not  account  for  radiation  losses.  The  scattering  parameters  determined  by  Touchstone 
satisfy  the  relationship  S,2  +  =  1.  The  circuit  simulator  results  provide  a  decent 

benchmark,  and  also  show  the  value  of  using  a  full-wave  solution.  There  is  some 
discrepancy  in  the  phase  of  the  circuit  simulator  compared  with  the  FDTD  results. 
However,  this  is  most  likely  because  there  is  radiation  from  the  line  which  is  not 
accounted  for  by  Touchstone.  There  is  no  difference  between  the  scattering  paramteters 
found  using  one  FDTD  run  and  two  FDTD  runs.  Using  the  Prony  technique  for 
discontinuity  problems  results  in  a  savings  of  50%  in  simulation  time  without  loss  of 
accuracy.  Moreover,  this  method  can  be  used  for  accurate  calculation  of  the  dispersive 
characteristics  of  microstrip  lines,  as  demonstrated  in  Section  4.2. 


Figure  4.7  Magnitude  of  the  scattering  parameters  of  the  right-angle  bend  as  a  function 
of  frequency. 


Frequency  (GHz) 


Figure  4.8  Phase  of  Sn  of  the  right-angle  bend  as  a  function  of  frequency. 


Figure  4.9  Phase  of  S2,  of  the  right-angle  bend  as  a  function  of  frequency 


4.5  Microstrip  Line  with  Mitered  Right- Angle  Bend 

One  limitation  of  the  orthogonal  FDTD  method  is  that  since  the  computational 
volume  is  discretized  with  cells  shaped  like  bricks,  geometries  containing  features  with 
curves  or  angles  that  are  not  90°  cannot  be  modeled  precisely.  The  usual  practice  is  to 
use  a  staircase  approximation  to  the  geometry.  In  terms  of  radiation  and  scattering 
problems,  staircasing  usually  gives  good  results  since  the  current  distribution  in  the  near¬ 
field  need  not  be  computed  exactly  to  predict  the  far-field  accurately.  However,  when 
this  technique  is  applied  to  guided-wave  or  circuit  problems,  staircasing  does  have  a 
considerable  effect  on  the  results,  as  will  be  demonstrated. 

A  solution  to  the  problem  of  accurately  modeling  geometries  that  are  non- 
Cartesian  is  to  use  the  general  curvilinear  FDTD  method.  This  method,  however,  is 
expensive  in  terms  of  the  time  required  for  the  update  equations  and  the  memory 
required,  but  does  produce  very  accurate  results.  An  alternative  to  the  full  curvilinear 
approach  is  to  write  a  modified  update  equation  which  can  be  applied  locally  to  the  cells 
which  contain  non-Cartesian  features.  In  order  to  demonstrate  this  approach,  consider  the 
infinitely  thin  perfect  electric  conductor  used  to  model  microstrip  lines.  Modeling  lines 
that  run  at  an  angle  with  respect  to  the  underlying  Cartesian  coordiante  system,  as  well  as 
lines  with  mitered  corners,  are  examples  of  problems  well-suited  to  the  modified  update 
equation. 

Another  common  microstrip  discontinuity  is  the  mitered  right-angle  bend,  which 
is  shown  in  Fig.  4.10.  The  45°  miter  is  usually  approximated  in  a  staircase  fashion  as 
shown.  A  better  approach  is  to  model  the  metal  with  a  triangular  cell.  Starting  with 
Faraday’s  Law  in  integral  form, 


(4.5) 


and  then  writing  the  discretized  version  of  the  equation  for  Hx  result  in 


76 


Hr'%,k)=H;-'%,k)+-^[E;(j,k)-E;u,k-i)\ 

—[^(J.khE'AJ-U)]  (4.6) 

Now,  assume  that  the  magnetic  field  H'  is  positioned  at  the  center  of  the  triangle  as 
shown  in  Fig.  4.10.  Following  the  integration  path  as  shown  in  the  figure  gives  the 
following  update  equation. 

H?'l2(j,k)=H;-'*(j,k)  +  -^{E;(j,k)Ay+  E;0-U)Az]  (4.7) 

Equation  (4.7)  was  implemented  to  model  infinitely  thin  triangular  pec  sections. 


Mitered  Right- Angle  Bend 


Figure  4.10  Mitered  right-angle  bend  with  detailed  FDTD  cell. 


The  mitered  right-angle  bend  was  simulated  with  the  staircase  approximation  and 
the  triangle  cell  approximation.  The  mesh,  excitation,  and  reference  planes  are  the  same 
as  shown  in  Fig.  4.6.  To  show  the  effect  of  the  miter  on  the  bend,  Fig.  4.1 1  shows  the 
voltage  at  reference  plane  1.  Here  the  triangle  approximation,  the  staircase 
approximation,  and  the  square  right-angle  bend  cases  are  compared.  The  miter  clearly 
reduces  the  reflection,  and  there  is  a  slight  difference  between  the  two  miter 
approximations.  Figures  4.12,  4.13,  and  4.14  compare  the  scattering  parameters  found 
via  FDTD  for  the  mitered  bend  with  those  calculated  using  Touchstone.  While  there  is 
not  much  difference  between  the  values  of  S2i  computed  using  the  two  approximations, 
there  is  a  significant  difference  when  comparing  the  Sn  values.  Although  the  staircase 
approximation  is  very  fine,  it  results  in  a  value  of  S„  that  is  larger  than  the  EEsof  result. 
On  the  other  hand,  the  triangular  cell  approximation  gives  extremely  good  agreement 
with  EEsof  up  to  about  9  GHz. 
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Figure  4.11  Voltage  at  reference  plane  1  as  a  function  of  time. 
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Figure  4.12  Magnitude  of  scattering  parameters  for  the  mitered  right-angle  bend  as  a 
function  of  frequency. 
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Figure  4.13  Phase  of  Su  for  the  mitered  right-angle  bend  as  a  function  of  frequency. 
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Figure  4.14  Phase  of  S2]  for  the  mitered  right-angle  bend  as  a  function  of  frequency. 

4.6  Skewed  Stripline 

The  effect  of  the  staircase  on  the  fields  propagating  along  a  microstrip  line  can  be 
determined  by  analyzing  a  stripline  that  runs  at  an  angle  of  45°  with  respect  to  the 
underlying  Cartesian  grid.  The  geometry  and  approximations  are  shown  in  Fig.  4.15, 
where  W  =  2.07  mm,  H  =  5.08  mm,  and  the  relative  dielectric  constant  is  5.  The  stripline 
structure  was  excited  with  a  Blackman-Harris  window  function  with  a  3  dB  cutoff 
frequency  of  8.65  GHz,  and  the  voltage  was  monitored  at  several  locations  along  the  line. 
An  expression  for  the  effective  dielectric  constant  in  terms  of  the  Fourier  transform  of  the 
voltage  at  two  locations,  separated  by  a  distance  L  on  the  line,  is  given  in  (4.8). 


where  c0  is  the  speed  of  light  in  free  space,  and  V(d)  is  the  Fourier  transform  of  the 


voltage  at  position  d.  Since  the  structure  being  analyzed  is  stripline,  it  is  known  that  the 
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Figure  4.15  (a)  Cross-section  of  stripline  with  W  =  2.07  mm  and  H  =  5.08  mm.  (b)  Line 
at  angle  of  45°  with  respect  to  the  underlying  grid,  (c)  Staircase 
approximation  to  the  line,  (d)  Triangular  cell  approximation  to  the  line. 


value  of  ereff{(0 )  is  a  constant,  and,  in  this  particular  case,  the  constant  is  5. 
Equation  (4.8)  will  be  used  to  calculate  ereff  for  a  straight  stripline  and  for  the  45°  line 
with  the  staircase  and  triangular  cell  approximations.  The  deviation  of  ereff  from  the 
expected  value  will  gauge  the  effectiveness  of  the  approximations. 
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The  voltage  as  a  function  of  time  at  two  locations  7.7625  mm  apart  on  a  straight 
stripline  is  shown  in  Fig.  4.16  Note  that  there  is  virtually  no  dispersion  of  the  voltage 
pulse  as  it  propagates  along  the  line.  Figure  4.17  shows  the  voltage  as  a  function  of  time 
at  two  locations  10.98  mm  apart  on  the  angled  stripline  of  Fig.  4.15.  There  is  a 
considerable  amount  of  dispersion  in  both  approximations.  Since  the  voltage  values  are 
similar,  it  is  helpful  to  examine  the  difference  between  the  two.  The  voltage  from  the 
triangular  cell  approximation  is  subtracted  from  the  staircase  voltage  and  plotted  in 
Fig.  4.18.  This  clearly  shows  that  the  voltage  travels  slower  when  the  staircase 
approximation  is  used.  Exactly  how  much  slower  the  voltage  travels  can  be  ascertained 
by  examining  Fig.  4.19  which  is  a  plot  of  the  relative  effective  dielectric  constant  as  a 
function  of  frequency.  The  staircase  approximation  yields  an  error  of  about  9%,  whereas 
the  triangular  cell  approximation  gives  an  error  of  at  most  3%.  This  is  significant,  since 
doubling  the  mesh  density  using  the  staircase  approximation  still  gives  an  error  of  about 
5%  [4.8].  Using  the  triangular  cell  approximation,  the  accuracy  is  improved  by  a  factor 
of  three  without  increasing  the  size  of  the  problem  or  slowing  the  update  equations. 
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Figure  4.16  Voltage  waveforms  as  a  function  of  time  at  two  positions  along  a  straight 
stripline.  The  monitoring  positions  are  7.7625  mm  apart. 


Figure  4.17  Voltage  waveforms  as  a  function  of  time  at  two  positions  along  an  angled 
stripline.  The  monitoring  positions  are  10.98  mm  apart. 
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Figure  4.18  Difference  between  the  staircase  and  triangular  cell  approximations  in  the 
voltage  waveforms  of  Fig.  4.17.  The  approximations  give  rise  to  different 
propagation  constants. 
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Figure  4.19  Effective  dielectric  constant  as  a  function  of  frequency.  The  triangular  cell 
approximation  is  three  times  more  accurate  than  the  staircase  approximation. 


4.7  Straight  Stripline 

In  the  previous  section,  a  skewed  stripline  was  analyzed.  The  effective  dielectric 
constant  was  calculated  as  a  function  of  frequency  from  the  time  domain  voltage.  In 
order  to  limit  the  source  of  error  to  the  staircase  or  triangular  cell  approximations,  the 
simulations  were  stopped  before  reflections  from  the  ABCs  could  contaminate  the  time 
domain  data.  In  this  section,  the  straight  stripline  is  again  considered.  However,  in  this 
section,  the  reflections  from  the  ABCs  terminating  the  computational  domain  will  have 
an  effect  on  the  calculation  of  the  effective  dielectric  constant.  Two  different  ABCs  were 
used  on  the  end  walls  of  the  stripline,  the  first-order  Mur  and  the  DBC.  The  parameters 
for  the  DBC  were  chosen  as  e  =  eh  =5.0  and  8  =  0.01.  The  magnitude  of  the 
reflection  coefficient  due  to  the  ABCs  is  shown  in  Fig.  4.20.  Both  ABCs  work  very  well 
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absorbing  the  TEM  wave  propagating  along  the  stripline.  The  reflection  is  less  than 
0.02%  for  the  DBC  in  the  0.5  GHz  to  4.5  GHz  frequency  range. 

Figure  4.21  shows  the  effective  dielectric  constant  as  calculated  using  (4.8).  The 
reflection  from  the  first-order  Mur  causes  oscillations  in  the  effective  dielectric  constant 
results.  The  oscillatory  behavior  is  not  present  in  the  effective  dielectric  constant  curve 
calculated  using  the  voltage  found  with  the  DBC.  The  Mur  and  DBC  results  are 
compared  to  a  case  in  which  no  reflection  is  present  in  the  time  domain  data.  When  no 
reflection  is  present,  the  time  domain  voltage  signature  is  four  times  shorter  than  that  with 
the  reflection  present.  This  explains  the  difference  in  the  frequency  resolution  between 
the  Mur,  DBC,  and  straight  line  cases.  Figure  4.22  shows  the  effective  dielectric  constant 
calculated  using  the  Prony  technique  of  Section  4.3.  The  Prony  technique  eliminates  all 
oscillatory  behavior  from  the  curves.  The  results  are  again  compared  with  those  for  the 
straight  line  case,  which  does  not  have  a  reflection  present  in  the  time  domain  data. 


Frequency  (GHz) 


Figure  4.20  Magnitude  of  reflection  coefficients  of  the  absorbing  boundary  conditions 
applied  to  the  straight  stripline  problem.  Both  the  first-order  Mur  and  the 
DBC  perform  extremely  well  for  this  case. 
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Frequency  (GHz) 


Figure  4.21  Effective  dielectric  constant  for  the  straight  stripline.  The  effect  of  the 
absorbing  boundaries  is  shown.  The  "straight  line"  curve  was  obtained  from 
a  time  signature  with  no  reflection  present. 


Figure  4.22  Effective  dielectric  constant  for  the  straight  stripline.  The  Prony  technique 
eliminates  the  effect  of  the  absorbing  boundary.  The  "straight  line"  curve 
was  obtained  from  a  time  signature  with  no  reflection  present. 


4.8  Coplanar  Waveguide  to  Microstrip  Line 

Designing  and  constructing  practical  microwave  and  millimeter-wave  circuits 
often  require  an  accurate  characterization  of  discontinuities  that  cannot  be  adequately 
modeled  using  circuit  simulators.  This  is  especially  true  at  higher  frequencies  where  a 
full-wave  solution  to  Maxwell's  equations  is  needed.  In  certain  circumstances,  when  two 
different  types  of  waveguiding  structures  are  connected,  there  are  simply  no  models  of 
these  types  of  transitions  in  circuit  simulators.  Some  examples  of  transitions  typically  not 
found  in  circuit  simulators  are  the  transitions  of  a  microstrip  line  to  a  grounded  coplanar 
waveguide,  a  slotline  to  a  finline  waveguide,  and  a  microstrip  line  to  a  waveguide.  The 
analysis  of  the  transition  from  a  coplanar  waveguide  to  a  microstrip  line  is  presented  in 
this  section. 

Both  quasi-static  and  full-wave  solution  methods  have  been  used  to  characterize  a 
coplanar  waveguide  [4.9].  A  coplanar  waveguide  is  preferred  over  a  microstrip  line  for 
several  reasons  [4.10].  Since  the  electric  field  is  confined  to  the  narrow  slot  region,  a 
relatively  small  amount  of  the  field  is  subjected  to  substrate  losses  as  compared  with  the 
microstrip  line.  Additionally,  having  the  ground  and  signal  lines  in  the  same  plane  is 
desirable  when  including  devices  in  the  circuit.  When  the  ground  and  signal  lines  are  not 
on  the  same  plane,  a  hole  through  the  substrate  is  required  to  reach  electric  ground;  this 
situation  is  impractical.  Moreover,  coupling  effects  between  adjacent  lines  are  very  weak 
since  a  ground  plane  separates  them.  When  an  additional  ground  plane  is  present,  as  is 
the  case  with  a  grounded  coplanar  waveguide,  the  gap  width  should  be  small  compared  to 
the  substrate  thickness  in  order  to  minimize  the  effect  of  the  microstrip  mode. 

The  top  view  of  the  transition  from  a  grounded  coplanar  waveguide  to  a 
microstrip  line  is  shown  in  Fig.  4.23.  The  thickness  of  the  dielectric  separating  the 
ground  planes  is  254  pm  and  the  dielectric  constant  is  2.2.  The  thickness  of  the  top 
ground  plane  and  microstrips  is  16.9333  pm.  The  width  of  the  conducting  strip  of  the 
coplanar  waveguide,  Wl,  is  442  pm,  and  the  width  of  the  slot,  SI,  is  54.5  pm.  The 
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transition  from  the  coplanar  waveguide  to  the  microstrip  line  is  1.1  mm  long,  and  the 
taper  angle  of  the  top  ground  plane  is  45°.  The  distance  from  the  top  ground  plane  to  the 
microstrip  line,  S2,  is  930  pm,  the  taper  angle  of  the  center  conductor  is  1 1.5°,  and  the 
width  of  the  microstrip  line  is  889  pm. 

The  mesh  used  to  describe  this  problem  is  very  dense  because  of  the  finite 
thickness  of  the  conducting  strips  and  the  need  to  accurately  model  the  slot  in  the 
coplanar  waveguide.  The  mesh  was  discretized  with  Ax  =  50  pm,  Ay  =  8.5  pm,  and 
A z  =  16.933  pm.  A  nonuniform  mesh  was  used  in  the  y-direction  in  order  to  have  at 
least  3*W2  of  space  between  the  microstrip  line  and  the  end  of  the  mesh.  The  largest  Ay 
was  17  pm.  First-order  Mur  ABCs  were  used  on  the  y-walls  and  the  top  wall.  The 
bottom  wall  was  pec,  and  the  z-walls  were  terminated  with  DBCs  with  £  =  1.403, 

er  =  2.1,  and  5=0.01. 

r2  ’ 


Figure  4.23  Top  view  of  the  grounded  coplanar  waveguide  to  microstrip  line  transition. 

W 1=442  pm,  W2=889  pm,  Sl=54.5  pm,  and  S2=930  pm,  and  the  transition 
is  1.1  mm  long. 
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The  mesh  had  40  cells  in  the  x-direction,  300  cells  in  the  y-direction,  and  152 
cells  in  the  z-direction.  The  discontinuity  was  22  cells  long  in  the  z-direction  and 
centered  in  the  mesh.  The  electric  field  in  the  slot  of  the  coplanar  waveguide  was  excited 
with  a  Blackman-Harris  window  function  with  a  10%  cutoff  frequency  of  100  GHz.  The 
voltage  in  the  slot  and  the  current  around  the  center  conductor  were  monitored  at  several 
locations  along  the  coplanar  waveguide.  The  voltage  under  the  center  conductor  was 
monitored  in  the  transition  region,  and  on  the  microstrip  side  of  the  problem,  the  voltage 
and  current  were  monitored  at  several  locations.  Time  domain  plots  of  the  voltage  are 
shown  in  Figs.  4.24  -  4.26.  The  numbers  in  the  plot  legends  refer  to  the  cell  monitoring 
location  in  the  z-direction.  The  voltage  plots  in  these  figures  show  that  the  voltage 
undergoes  considerable  reflection  because  of  the  transition  as  well  as  dispersion  due  to 
the  air-dielectric  interface. 

Using  the  Prony  technique  to  extract  forward  and  backward  traveling  voltage  and 
current  waves,  the  impedance  of  the  lines  was  calculated  and  is  shown  in  Fig.  4.27.  The 
coplanar  waveguide  impedance  is  a  few  ohms  lower  than  the  desired  50  £2  at  77  GHz. 


Figure  4.24  Voltage  in  the  slot  of  the  grounded  coplanar  waveguide  as  a  function  of 
time.  The  voltage  is  shown  at  three  successive  monitoring  locations  in  the 
z-direction. 
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Figure  4.25  Voltage  under  the  conductor  in  the  transition  region  as  a  function  of  time. 

The  voltage  is  shown  at  four  successive  monitoring  locations  in  the 
z-direction. 
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Figure  4.26  Microstrip  line  voltage  as  a  function  of  time.  The  voltage  is  shown  at  three 
successive  monitoring  locations  in  the  z-direction. 
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Figure  4.27  Magnitude  of  the  impedance  of  the  grounded  coplanar  waveguide  and  the 
microstrip  line  as  a  function  of  frequency. 


Since  the  impedances  on  either  side  of  the  transition  discontinuity  are  not  the 
same,  normalized  scattering  parameters  must  be  used.  The  expression  for  S21  is  given  in 


(4.9) 


S21(f)  = 


[MUMP 

W)  /,(/) 


(4.9) 


where  I2(f)  is  the  Fourier  transform  of  the  time  domain  current  on  the  microstrip  side  of 
the  discontinuity,  Z2(f)  is  the  impedance  of  the  microstrip  line,  and  the  subscript  "1" 
refers  to  the  grounded  coplanar  waveguide.  Running  the  simulation  again  with  the 
excitation  on  the  microstrip  line  side  of  the  discontinuity  enables  the  calculation  of  all 
four  scattering  parameters.  The  magnitudes  of  the  scattering  parameters  of  the  coplanar 
waveguide  to  microstrip  line  transition  are  shown  in  Fig.  4.28.  The  phase  is  shown  in 
Fig.  4.29.  It  should  be  noted  that  S21  and  S;2  were  slightly  different,  so  an  average  of  the 
two  values  was  taken,  thereby  enforcing  reciprocity. 


Figure  4.28  Magnitude  of  the  scattering  parameters  as  a  function  of  frequency  for  the 
grounded  coplanar  waveguide  to  microstrip  line  transition. 


Figure  4.29  Phase  of  the  scattering  parameters  as  a  function  of  frequency  for  the 
grounded  coplanar  waveguide  to  microstrip  line  transition. 
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A  final  interesting  quantity  to  examine  is  the  normalized  power  loss,  which  is 
given  by 

Power  Loss  =  1  - 15, ,  |2  - 1^,  |2  (4.1 0) 

when  the  source  is  on  the  port-one  side  of  the  discontinuity.  The  normalized  power  loss 
is  shown  as  a  function  of  frequency  in  Fig.  4.30.  The  PI  in  the  legend  refers  to  the  source 
on  the  coplanar  waveguide  side,  and  the  P2  refers  to  the  source  on  the  microstrip  line 
side. 


Figure  4.30Normalized  power  loss  as  a  function  of  frequency  for  the  coplanar 
waveguide  to  microstrip  line  transition.  PI  refers  to  the  source  on  the 
coplanar  waveguide  side,  and  P2  refers  to  the  source  on  the  microstrip  line 
side. 


4.9  Conclusions 

A  variety  of  microstrip  line  and  stripline  structures  were  characterized  in  this 
chapter  using  the  FDTD  method.  A  technique  based  on  Prony's  method  and  transmission 
line  equations  was  introduced,  and  demonstrated  to  characterize  lines  by  extracting 
reflections  from  ABCs  and  to  enable  one  to  calculate  scattering  parameters  of  a  two-port 
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system  with  a  single  FDTD  simulation.  Because  of  the  expense  of  FDTD  simulations  in 
terms  of  memory  and  computation  time,  reducing  the  number  and  length  of  simulations  is 
a  worthwhile  goal.  A  special  update  equation  based  on  a  triangular  path  of  integration 
was  introduced  and  shown  to  be  a  significant  improvement  over  the  conventional 
staircase  method  for  approximating  planar  structures  that  do  not  directly  correspond  to 
the  underlying  Cartesian  grid.  Finally,  the  transition  from  grounded  coplanar  waveguide 
to  microstrip  line  was  investigated.  Both  the  nonuniformity  of  the  mesh  and  the  Prony 
technique  for  scattering  parameter  extraction  were  important  in  this  example.  The 
grounded  coplanar  waveguide  design  presented  very  nearly  achieves  a  match  at  77  GHz. 
However,  there  is  still  significant  power  loss.  While  it  is  believed  that  a  certain  amount 
of  power  loss  is  inevitable  due  to  the  geometry  of  the  discontinuity,  perhaps  the  power 
loss  could  be  abated  through  the  use  of  a  tuning  stub  on  the  coplanar  waveguide. 
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CHAPTER  5 

CHARACTERIZATION  OF  COUPLED  MICROSTRIP  LINES  AND  COUPLED 
STRIPLINES  USING  NONUNIFORM  ORTHOGONAL  FDTD 

5.1  Introduction 

Coupled  microstrip  lines  and  coupled  striplines  are  used  in  the  design  of 
directional  couplers,  filters,  impedance  matching  networks,  baluns,  and  delay  lines  [5.1]. 
Coupled  lines  are  also  important  in  digital  circuits  where  undesired  crosstalk  can  occur 
when  signal  lines  are  brought  into  close  proximity  of  each  other  [5.2].  There  has  been 
extensive  work  published  in  the  area  of  coupled  line  analysis  and  design.  A  selected 
sample  of  some  of  this  work  is  listed  in  the  references,  [5.3]  -  [5.1 1].  Various  numerical 
methods  have  been  employed  in  the  study  of  coupled  lines,  including  the  even  and  odd 
mode  method,  the  coupled  mode  formulation,  the  Green's  function  method,  the 
variational  method,  Fourier- series  expansion  method,  and  the  conformal  mapping 
technique  [5.4].  In  [5.7]  the  spectral  Galerkin  procedure  was  used.  Time  domain  models 
were  developed  in  [5.6]  and  [5.8]  -  [5.10]  utilizing  quasi-TEM  matrices.  These  methods 
require  the  solution  of  an  eigenvector  problem.  The  asymmetric  coupled  line  case  is 
defined  and  treated  with  the  quasi-TEM  approach  in  [5.10]. 

The  purpose  of  this  chapter  is  to  use  the  full-wave  solution  provided  by  the  FDTD 
method  to  characterize  coupled  microstrip  lines  and  coupled  striplines.  First,  the 
symmetric  problem  is  considered  and  an  even  and  odd  mode  theory  is  used  in  addition  to 
the  methods  of  Chapter  4  to  characterize  the  lines.  The  methods  of  Chapter  4  are 
extended  to  handle  the  asymmetric  case.  This  new  method  is  presented,  followed  by 
numerical  examples  analyzing  coupled  symmetric  and  asymmetric  lines.  Where  possible, 
FDTD  results  are  compared  with  circuit  simulator  data.  Finally,  a  multidielectric 
asymmetric  coupled  stripline  problem  is  analyzed. 


5.2  Symmetric  Coupled  Microstrip  Lines 

In  this  section,  the  symmetric  coupled  microstrip  lines  of  Fig.  5.1  are  analyzed 
using  FDTD  and  the  Prony  technique  described  in  Section  4.3.  The  lines  are  parallel  and 
run  in  the  z-direction.  The  width  of  each  line  is  300  (im  and  the  separation  between  the 
lines  is  also  300  pm.  The  substrate  below  the  infinitely  thin  lines  is  250  pm  thick  and  has 
a  dielectric  constant  of  4.5. 

The  general  N-line  system  supports  N  different  modes.  For  this  case,  there  are 
two  lines,  so  there  are  two  modes.  Because  of  the  symmetry  of  the  structure,  it  is 
straightforward  to  determine  the  field  distributions  of  the  two  modes,  which  can  be 
referred  to  as  the  even  and  odd  modes.  For  the  even  mode,  the  lines  have  the  same 
voltage,  whereas  for  the  odd  mode,  the  lines  have  voltages  which  are  equal  in  magnitude 
but  180  degrees  out  of  phase  [5.1].  One  approach  for  determining  the  effective  dielectric 
constant  and  impedance  of  each  mode  is  to  model  the  coupled  lines  twice,  once  with  an 
even  mode  excitation,  and  once  with  an  odd  mode  excitation.  Using  voltage  and  current 
data,  the  effective  dielectric  constant  and  impedance  are  found  as  detailed  in  Chapter  4. 
Another  approach  is  to  use  a  symmetry  plane  between  the  lines.  Again  the  problem  is 
solved  twice,  once  using  a  plane  of  perfect  electric  conductor  for  the  odd  mode,  and  once 
using  a  plane  of  perfect  magnetic  conductor  for  the  even  mode  [5.1 1]. 


Figure  5.1  Geometry  of  coupled  microstrip  lines.  For  the  symmetric  case  in  this 
section,  S  =  W1  =  W2  =  300  pm.  The  substrate  thickness,  H,  is  250  pm. 


The  method  used  in  this  work  requires  only  one  FDTD  simulation.  The  lines  are 
excited  with  an  arbitrary  voltage  distribution,  which  by  definition  is  a  linear  superposition 
of  the  even  and  odd  mode  field  distributions  [5.1].  For  the  example  presented  in  this 
section,  line  1  was  excited  and  line  2  was  not  excited.  Figure  5.2  shows  the  line  current 
on  line  1  at  two  positions  on  the  line  as  a  function  of  time.  Position  1  is  located  0.25  cm 
from  the  source  and  position  2  is  1.125  cm  from  the  source.  The  amplitude  of  the  current 
on  line  1  decays  as  energy  is  coupled  into  line  2  as  the  voltage  and  current  travel  down 
the  line.  This  phenomenon  is  clearly  displayed  in  Fig.  5.3,  which  shows  the  line  current 
on  line  2  as  a  function  of  time.  In  Fig.  5.3  the  amplitude  of  the  current  increases  as  the 
distance  from  the  source  increases. 

The  even  and  odd  mode  currents  are  given  in  (5.1)  and  (5.2). 

Ieven={li+I2)/2.0  (5.1) 

I  odd  ~{h  ~/2)/2-0  (5.2) 

where  7,  and  /2  are  the  currents  on  line  1  and  line  2,  respectively.  Equations  (5.1)  and 
(5.2)  are  valid  in  both  the  time  domain  and  the  frequency  domain.  The  even  mode 
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Figure  5.2  Line  current  on  line  1  as  a  function  of  time.  Position  1  and  position 
2  are  0.875  cm  apart. 


current  was  determined  using  (5.1),  and  is  plotted  as  a  function  of  time  in  Fig.  5.4. 
Equation  (5.2)  was  used  to  calculate  the  odd  mode  current,  which  is  shown  in  Fig.  5.5. 
Comparing  Fig.  5.4  and  Fig.  5.5,  it  is  clear  that  the  odd  mode  current  propagates  faster 
than  the  even  mode  current. 


Time  Step 


Figure  5.3  Line  current  on  line  2  as  a  function  of  time.  Position  1  and  position  2  are 
0.875  cm  apart. 
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Figure  5.4  Even  mode  current  as  a  function  of  time.  Position  1  and  position  2  are 
0.875  cm  apart. 
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Figure  5.5  Odd  mode  current  as  a  function  of  time.  Position  1  and  position  2  are 
0.875  cm  apart. 

The  effective  dielectric  constant  and  impedance  are  found  using  the  procedure  in 
Section  4.3,  which  takes  into  account  the  presence  of  reflections  in  the  time  domain  data 
due  to  the  ABCs.  The  even  and  odd  mode  effective  dielectric  constants  are  shown  in 
Fig.  5.6.  FDTD  results  are  compared  with  results  from  the  circuit  simulator  LINECALC 
[5.12].  The  agreement  for  the  odd  mode  is  excellent.  For  the  even  mode,  however,  there 
is  about  a  4%  difference  between  the  FDTD  and  LINECALC  values  at  the  low  frequency 
end  of  the  graph.  The  difference  is  less  than  3%  at  a  frequency  of  40  GHz.  The 
impedance  is  plotted  as  a  function  of  frequency  in  Fig.  5.7.  Again,  the  FDTD  results  are 
compared  with  the  LINECALC  results.  The  agreement  between  the  FDTD  values  and 
the  LINECALC  values  is  excellent.  The  impedance  of  a  single  isolated  line  as  calculated 
using  LINECALC  is  provided  as  a  reference.  At  low  frequencies,  the  even  mode 
impedance  is  about  12%  larger  than  the  single  line  impedance.  At  high  frequencies,  the 
even  mode  impedance  is  about  15%  larger  than  the  single  line  impedance.  The  odd  mode 
impedance  is  about  14%  smaller  than  the  single  line  impedance. 


Figure  5.6  Effective  dielectric  constant  as  a  function  of  frequency  for  the  symmetric 
coupled  microstrip  lines  in  Fig.  5.1. 


Figure  5.7  Even  and  odd  mode  impedances  as  functions  of  frequency  for  the  symmetric 
coupled  microstrip  lines  in  Fig.  5 . 1 . 
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5.3  Prony's  Method  for  Mode  Extraction 

A  new  approach  for  mode  extraction  of  coupled  lines  is  presented  in  this  section. 
Neglecting  the  reflections  from  the  ABCs,  the  voltage  on  either  line  of  a  coupled  two-line 
system  can  be  written  in  terms  of  the  modal  voltages  as 

V(z)  =  (5.3) 

where 

r,=«,+;A  (5.4) 

and  A,  is  the  complex  amplitude  of  the  voltage  mode.  Similarly,  an  expression  for  the 
current  is  given  in  (5.5) 


C(z)  =  Bler't  +  Bjer*  (5.5) 

where  y,  is  defined  in  (5.4).  There  are  six  unknowns  in  (5.3)  and  (5.5),  namely  Ax,  A1, 
Bv  B2,  and  y2.  Once  these  unknowns  are  determined,  it  is  straightforward  to 
calculate  the  modal  effective  dielectric  constants  and  modal  impedances. 

The  unknowns  in  (5.3)  and  (5.5)  can  be  determined  using  Prony's  method  [5.13]. 
The  following  derivation  will  be  based  on  the  expression  for  the  voltage  in  Equation 
(5.3).  The  voltage  is  monitored  at  four  or  more  locations  in  the  z-direction.  These 
monitoring  points  should  be  equally  spaced.  The  four  time  domain  voltage  responses  are 
transformed  to  the  frequency  domain  and  written  as  V0,  Vx,  V2,  and  Vv  Equation  (5.3)  is 
rewritten  as  (5.6). 


The  m,  are  then  given  by 


where 


and 


V(z)  =  Aii^+A2i4 

(5.6) 

(5.7) 

A  _  YYi-VoVt 

(5.8) 

(5.9) 

and  the  two  values  of  w,  are  determined  by  the  ±  operator  in  (5.7). 
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When  applying  this  method,  it  is  important  that  reflections  from  the  ABCs  are  not 
present  in  the  time  domain  data.  Reflections  can  be  avoided  by  simply  making  the 
problem  long  in  the  z-direction  and  stopping  the  simulation  before  reflections  from  the 
ABCs  can  contaminate  the  fields  at  the  monitoring  locations.  This  will  reduce  the 
number  of  points  in  the  time  domain  which  leads  to  coarser  resolution  in  the  frequency 
domain.  For  many  applications,  this  is  not  a  problem,  since  the  frequency  domain  data 
are  smooth  and  do  not  require  fine  resolution. 

5.4  Coupled  Microstrip  Lines 

The  technique  presented  in  the  previous  section  will  be  demonstrated  in  this 
section  to  calculate  the  effective  dielectric  constant  and  impedance  for  the  two  modes  of 
the  coupled  microstrip  lines  shown  in  Fig.  5.1.  For  this  case,  W1  =  W2  =  300  |i.m,  and 
S  =  50  pm.  Since  the  lines  are  close  together,  it  is  expected  that  the  coupling  effects  will 
be  greater  than  in  the  example  in  Section  5.2. 

Line  1  was  excited  with  a  Blackman-Harris  window  function  with  a  3  dB  cutoff 
frequency  of  78.9  GHz  and  line  2  was  not  excited.  The  voltage  and  current  were 
monitored  along  the  coupled  lines,  transformed  to  the  frequency  domain,  and  used  to 
calculate  effective  dielectric  constants  and  impedances.  Figure  5.8  shows  the  two  modal 
effective  dielectric  constants  and  compares  them  with  the  even  and  odd  effective 
dielectric  constants  calculated  using  LINECALC.  The  agreement  between  the 
LINECALC  and  FDTD  values  is  excellent  for  the  odd  mode.  For  the  even  mode,  the 
difference  is  only  1.3%.  Even  and  odd  mode  impedances  are  shown  in  Fig.  5.9.  The 
agreement  for  the  even  mode  is  excellent  but  the  odd  mode  calculated  using  FDTD  is 
considerably  lower  than  the  LINECALC  impedance.  Increasing  the  mesh  density  in  the 
cross-section  would  improve  the  FDTD  results.  As  in  Fig.  5.7,  the  impedance  of  a  single 
isolated  line  is  shown  for  comparison  with  the  even  and  odd  mode  impedances.  As  the 
lines  are  moved  closer  together,  the  coupling  becomes  stronger,  and  there  is  greater 


deviation  of  the  even  and  odd  mode  impedances  from  the  characteristic  impedance  of  the 
isolated  line. 
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Figure  5.8  Effective  dielectric  constant  as  a  function  of  frequency  for  symmetric 
coupled  microstrip  lines  with  S  =  50  p.m. 
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Figure  5.9  Even  and  odd  mode  impedances  as  functions  of  frequency  for  symmetric 
coupled  microstrip  lines  with  S  =  50  |im. 


5.5  Asymmetric  Coupled  Microstrip  Lines 

Consider  the  geometry  shown  in  Fig.  5.1.  If  W1  *  W2,  the  method  presented  in 
Section  5.2  for  calculating  the  modal  effective  dielectric  constants  and  impedances  cannot 
be  used.  However,  the  method  presented  in  Section  5.3  can  be  used.  Numerical  results 
will  be  presented  for  asymmetric  coupled  lines  with  W1  =  300  pm,  W2  =  200  pm,  and 
S  =  50  pm.  Line  1  was  excited  and  energy  coupled  into  line  2.  The  voltage  and  current 
were  monitored  at  several  equally  spaced  sampling  locations  along  the  lines.  Effective 
dielectric  constants  and  impedances  were  then  calculated,  as  described  in  the  previous 
section. 

Figure  5.10  shows  the  effective  dielectric  constant  as  a  function  of  frequency. 
Values  from  the  FDTD  simulation  are  compared  to  values  from  LINECALC.  Since 
LINECALC  does  not  solve  the  asymmetric  case,  two  sets  of  LINECALC  data  are  shown. 
One  set  is  for  symmetric  lines  with  width,  W  =  300  pm,  and  separation,  S  =  50  pm. 
These  data  are  labeled  "Even  300"  and  "Odd  300"  in  the  graph  legend  in  Fig.  5.10.  The 
"Even  200"  and  "Odd  200"  labels  refer  to  LINECALC  results  for  symmetric  lines  with 
width,  W  =  200  pm,  and  separation,  S  =  50  pm.  As  shown  in  Fig.  5.10,  the  values  of  the 
effective  dielectric  constant  for  mode  1  from  the  FDTD  simulation  are  between  those 
from  the  LINECALC  calculations.  This  is  the  expected  result.  If  it  is  assumed  that 
FDTD  values  for  mode  2  are  slightly  larger  than  they  should  be,  as  was  the  case  in 
Fig.  5.6  and  Fig.  5.8,  the  FDTD  data  are  consistent  with  the  results  from  LINECALC. 

Modal  impedance  values  are  displayed  in  Figs.  5.1 1  and  5.12  for  line  1  and  line  2, 
respectively.  The  FDTD  values  are  very  reasonable,  with  the  exception  of  mode  1  on  line 
1,  which  is  most  likely  smaller  than  it  should  be.  However,  the  overall  results  suggest 
that  the  approach  presented  in  Section  5.3  is  a  good  method  for  characterizing 
asymmetric  coupled  microstrip  lines. 
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Figure  5.10  Even  and  odd  mode  effective  dielectric  constants  as  functions  of  frequency 
for  asymmetric  coupled  microstrip  lines  with  W1  =  300  (im,  W2  =  200  pm, 
and  S  =  50  pm.  The  FDTD  results  are  compared  with  circuit  simulator 
results  for  symmetric  coupled  lines  with  a  separation  of  50  pm. 
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Figure  5.11  Impedance  as  a  function  of  frequency  for  asymmetric  coupled  lines.  The 
modal  impedances  calculated  on  line  1  (W1  =  300  pm)  are  compared  with 
circuit  simulator  results  for  symmetric  coupled  microstrip  lines  with  width 
300  pm.  The  single  line  result  shows  the  impedance  of  an  isolated  line  with 
width  300  pm. 
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Figure  5.12  Impedance  as  a  function  of  frequency  for  asymmetric  coupled  microstrip 
lines.  The  modal  impedances  calculated  on  line  2  are  compared  with  circuit 
simulator  results  for  the  symmetric  case  with  W  =  200  pm.  The  single  line 
result  shows  the  impedance  of  an  isolated  line  with  width  200  pm 


5.6  Asymmetric  Coupled  Striplines 

The  previous  sections  in  this  chapter  dealt  with  the  analysis  of  coupled  microstrip 
lines,  and  FDTD  results  were  compared  with  results  from  LINECALC.  In  this  section, 
the  effective  dielectric  constant  of  asymmetric  coupled  striplines  will  be  determined. 
LINECALC  results  are  not  available  for  this  structure  because  of  the  asymmetry  and 
variation  of  the  dielectric.  The  geometry  is  shown  in  Fig.  5.13,  where  W1  =  300  pm, 
W2  =  200  pm,  S  =  50  pm,  HI  =  H2  =  250  pm,  and  the  top  and  bottom  layers  of  the 
dielectric  have  dielectric  constants  of  2.2  and  4.5,  respectively.  In  order  to  extract  the 
effective  dielectric  constant  for  the  asymmetric  coupled  striplines,  line  1  is  excited  with  a 
Blackman-Harris  window  function  with  a  3  dB  cutoff  frequency  of  45.6  GHz  and  line  2 
is  left  unexcited.  Figure  5.14  shows  the  line  current  on  line  1  as  a  function  of  time.  The 
amplitude  of  the  current  decays  as  energy  is  coupled  into  line  2.  This  is  shown  in 


Fig.  5.15,  where  the  amplitude  of  the  line  current  on  line  2  increases  as  the  monitoring 
location  moves  farther  away  from  the  source. 


Figure  5.13  Geometry  of  asymmetric  coupled  striplines.  W1  =  300  |im,  W2  =  200  |im, 
S  =  50  |lm,  and  HI  =  H2  =  250  |im. 
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Figure  5.14  Line  current  on  line  1  for  the  asymmetric  coupled  stripline  in  Fig.  5.13. 
Position  1  and  position  2  are  0.750  cm  apart. 
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Figure  5.15  Line  current  on  line  2  for  the  asymmetric  coupled  stripline  in  Fig.  5.13. 
Position  1  and  position  2  are  0.750  cm  apart. 


Figure  5.16  Effective  dielectric  constant  as  a  function  of  frequency  for  the  asymmetric 
coupled  stripline  in  Fig.  5.13.  The  modal  effective  dielectric  constants  are 
compared  with  the  effective  dielectric  constants  of  the  isolated  stripline 
embedded  in  the  mixed  dielectric.  W1  =  300  |im,  and  W2  =  200  jim. 
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The  effective  dielectric  constant  is  plotted  as  a  function  of  frequency  in  Fig.  5.16. 
The  two  modal  effective  dielectric  constants  are  compared  to  the  effective  dielectric 
constants  for  single  isolated  striplines  of  widths  300  pm  and  200  pm,  respectively.  Since 
LINECALC  results  are  not  available,  FDTD  results  for  the  isolated  lines  are  used  to 
compare  with  the  modal  results.  As  is  the  case  with  microstrip  lines,  the  modal  effective 
dielectric  constants  are  greater  than  and  less  than  the  effective  dielectric  constant  for  the 
single  line  case.  This  type  of  behavior  suggests  that  these  results  are  reasonable. 

5.7  Conclusions 

Coupled  microstrip  lines  and  coupled  striplines  were  analyzed  in  this  chapter. 
The  FDTD  method  was  used  in  conjunction  with  techniques  based  on  Prony's  method  to 
calculate  modal  effective  dielectric  constants  and  impedances  with  a  single  FDTD 
simulation.  Moreover,  the  new  approach  presented  was  not  dependent  upon  the 
symmetry  of  the  lines.  Numerical  results  for  symmetric  and  asymmetric  lines  were 
presented  and  compared  with  results  from  LINECALC.  An  asymmetric  coupled  stripline 
problem  with  mixed  dielectric  substrate  was  analyzed,  and  the  coupled  results  were 
compared  with  single  line  results  obtained  using  FDTD  since  problems  with  layered 
dielectrics  cannot  be  modeled  using  LINECALC.  The  results  from  the  method  presented 
were  shown  to  be  acceptable.  The  method  could  be  further  improved  to  account  for  the 
modal  reflections  due  to  the  ABCs.  This  would  require  determining  the  roots  of  a 
fourth-order  polynomial  rather  than  the  roots  of  a  quadratic.  Another  possible  avenue  of 
research  would  be  to  use  the  pencil  of  functions  method  [5.14]  to  set  up  the  eigenvector 
problem  for  the  general  coupled  N-line  system. 
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CHAPTER  6 

ANALYSIS  OF  MICROSTRIP  ANTENNAS  USING  FDTD 


6.1  Introduction 

The  FDTD  method  has  been  applied  to  analyze  a  variety  of  radiating  microstrip 
and  wire  structures  [6.1]-[6.10].  The  main  advantage  of  the  FDTD  method  over  moment 
method  techniques  for  complex  geometries  is  that  the  preprocessing  required  for  FDTD  is 
virtually  negligible  when  compared  to  the  preprocessing  required  to  derive  and 
implement  the  requisite  Green's  functions  for  the  method  of  moments.  In  this  chapter,  the 
near-field  to  far-field  transformation  is  discussed,  followed  by  the  analysis  of  a 
complicated  balun-fed  folded  dipole  antenna.  Next  the  benefits  of  using  the  nonuniform 
orthogonal  FDTD  method  are  demonstrated  through  the  analysis  of  microstrip  patch 
antennas.  The  results  obtained  are  compared  with  experimental  measurements. 

6.2  Near-field  to  Far-field  Transformation 

There  are  several  methods  for  obtaining  far-field  patterns  from  near-field  data  in 
an  FDTD  simulation  [6. 1 1]-[6. 13].  The  method  reported  in  [6.13]  provides  frequency 
results  over  a  wide  band,  but  requires  an  exorbitant  amount  of  memory  to  store  entire 
time  signatures  at  the  observation  locations.  In  this  work,  the  method  reported  in  [6.11] 
was  implemented.  It  is  a  single  frequency  method  that  is  not  expensive  in  terms  of 
memory  requirements  as  compared  to  the  method  in  [6.13].  A  brief  summary  of  the 
procedure  in  [6.11]  is  presented  here. 

To  determine  the  radiation  pattern  or  far-field  data  for  a  general  radiator  or 
scatterer,  one  typically  approaches  the  problem  by  solving  for  the  currents  on  the  surface 
of  the  radiator.  The  currents  are  then  transformed  to  the  far  field  using  a  Fourier 
transform.  However,  when  the  radiator  is  a  complex  geometry  containing  multiple  layers 
of  dielectric  and  metallization  in  both  the  vertical  and  horizontal  directions,  it  is 
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extremely  difficult  to  derive  an  accurate  Green's  function,  which  is  required  in  order  to 
solve  for  the  currents.  This  is  not  a  difficulty  for  the  FDTD  method,  which  solves  for  the 
fields  everywhere  within  the  computational  volume.  Therefore,  Huygens'  Equivalence 
Principle  can  be  applied  in  a  straightforward  manner.  The  radiator  is  enclosed  by  an 
equivalent  surface  as  shown  in  Fig.  6.1.  The  tangential  electric  and  magnetic  fields  on 
the  surface  are  used  to  calculate  the  equivalent  electric  and  magnetic  surface  currents. 
Electric  and  magnetic  potentials  are  calculated  from  the  currents  as  shown  in  Equation 


IffF  \ik°r'™lds’a 

Am  PI  m 


where  r'cos^  =  (jc,cos^  +  y/sin0)sin6  +  z/cos0.  Expressions  for  the  far  field  are 
determined  using  the  simple  free-space  Green's  function.  The  integration  is  carried  out 
numerically  by  a  straightforward  summation. 


Figure  6.1  Schematic  representation  of  radiating  source  surrounded  by  equivalent 
sources. 
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The  electric  field  componenents  Ee  and  are  then  given  by  Equations  (6.2)  and 


(6.3) 
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where  r\0  is  the  intrinsic  impedance  of  free  space.  It  should  be  noted  that  the  expressions 
in  (6.1)  -  (6.3)  are  written  in  the  frequency  domain.  Since  FDTD  is  a  time  domain 
method,  the  magnitude  and  phase  of  the  surface  currents  are  determined  by  assuming  that 
fields  have  settled  to  the  steady-state  solution.  The  maximum  and  minimum  values  of  the 
field  or  current  at  a  given  location  on  the  equivalent  surface  are  determined  by  searching 
the  data  and  the  slope  of  the  data  as  time  progresses.  However,  the  entire  time  signature 
is  not  stored.  Only  three  time  points  are  required  per  field  location. 


6.3  Baiun-Fed  Folded  Dipole 

The  objectives  of  the  work  presented  in  this  section  are  to  use  the  FDTD  method 
to  calculate  the  standing  wave  ratio  as  a  function  of  frequency,  and  to  calculate  far-field 
radiation  patterns  at  several  frequencies,  for  the  balun-fed  folded  dipole  designed  and 
tested  by  Proudfoot  [6.14].  Advantages  of  the  balun-fed  folded  dipole  were  reported  in 
[6.14]  and  will  be  summarized  here.  The  balun-fed  folded  dipole  combines  the  moderate 
bandwidth  of  the  folded  dipole  with  the  double  tuning  capabilites  of  the  balun  to  provide 
a  significant  increase  in  bandwidth  compared  to  that  of  conventional  microstrip  patches. 
The  high  impedance  of  the  folded  dipole  leads  to  a  narrow  microstrip  line  feed,  which  is 
important  from  the  fabrication  standpoint.  Low  impedances  lead  to  wide  lines,  which  in 
certain  instances  are  actually  wider  than  the  radiating  element.  A  final  advantage  of  the 
balun-fed  folded  dipole  is  that  its  endfire  properties  and  shielding  ground  plane  make  it 
attractive  for  use  in  arrays. 
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The  balun-fed  folded  dipole  antenna  shown  in  Figs.  6.2  and  6.3  is  a  fairly 
complicated  structure.  A  “tee”  shape  is  etched  away  from  copper  forming  the  conducting 
arms  of  the  folded  dipole.  The  length  of  the  folded  dipole  is  3.048  cm,  the  width  of  the 
top  arm  is  0.1016  cm,  the  width  of  the  bottom  arm  is  0.3810  cm,  and  the  etched  gap 
separating  the  arms  is  0.0508  cm.  The  dipole  is  excited  with  a  0.1151  cm  wide 
microstrip-line  balun  feed.  The  copper  of  the  folded  dipole  extends  towards  the  source,  is 
three  times  as  wide  as  the  microstrip  line,  and  serves  as  the  ground  for  the  feeding 
microstrip  line.  The  etched  slot  parallel  to  the  feed  line  has  the  same  width  as  the  feed 
line.  The  dipole  and  the  feeding  line  are  separated  by  a  dielectric  with  a  thickness  of 
0.16  cm  and  relative  permittivity  equal  to  4.4.  The  microstrip  line  is  impedance  matched 
at  the  orthogonal  ground  plane,  which  has  a  hole  in  it.  The  orthogonal  ground  plane  is 
40.64  cm  by  40.64  cm,  and  located  at  the  z  =  0  plane,  2.027  cm  from  the  center  of  the 
folded  dipole.  These  physical  dimensions  are  found  in  [6.14],  and  some  are  from 
measurements  of  an  element  at  the  Rome  Laboratory,  Hanscom  Air  Force  Base.  There  is 
some  uncertainty  associated  with  the  dimensions  determined  by  measurement  because  the 
measured  element  was  a  prototype. 

Examining  the  time  signature  of  the  voltage  gives  a  physical  picture  of  the  effects 
of  the  discontinuities.  Figure  6.4  shows  the  voltage  under  the  microstrip  line  at 
z  =  -1.0  cm  and  z  =  0.5  cm.  The  effects  of  the  hole  in  the  orthogonal  ground  plane  and 
the  antenna  are  shown  in  Fig.  6.4(a).  There  is  a  peak  in  the  voltage  response  due  to  the 
hole  at  time  step  400.  Later,  at  approximately  time  step  850,  the  reflection  caused  by  the 
antenna  is  seen.  Figure  6.4(b)  shows  that  there  is  a  difference  between  the  uniform  line 
case  and  the  antenna  case  before  the  discontinuity  of  the  antenna  is  encountered  because 
of  the  hole  in  the  orthogonal  ground  plane.  This  agrees  with  the  results  in  Fig.  6.4(a). 


Figure  6.2  Balun-fed  folded  dipole  antenna.  A  layer  of  dielectric  material  separates  the 
microstrip  line  feed  from  the  folded  dipole.  The  folded  dipole  is  an 
extension  of  the  ground  for  the  feed.  A  hole  in  the  orthogonal  ground  plane 
allows  the  microstrip  line  to  pass  through. 
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Figure  6.3  Cross-sections  of  the  balun-fed  folded  dipole  antenna,  (a)  The  x-y  plane, 
(b)  The  x-z  plane. 
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(a) 


(b) 

Figure  6.4  Voltage  along  the  microstrip  feed  line  as  a  function  of  time,  (a)  At 
z  =  -1.0  cm.  (b)  At  z  =  0.5  cm. 
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The  SWR  was  calculated  by  first  subtracting  the  time  signature  of  the  voltage  of 
the  uniform  line  from  the  time  signature  of  the  voltage  of  the  balun-fed  folded  dipole  with 
an  orthogonal  ground  plane  and  calling  this  the  reflected  voltage.  Next,  the  reflected 
voltage  and  the  incident  voltage  (the  voltage  of  the  uniform  line)  were  transformed  to  the 
frequency  domain  via  the  fast  Fourier  transform.  The  magnitude  of  the  SWR  was  then 
computed  with  the  following  formula: 


The  SWR  was  also  computed  using  the  magnitudes  of  the  incident  and  reflected  currents 
rather  than  the  magnitudes  of  the  incident  and  reflected  voltages  in  Equation  (6.4). 

In  Fig.  6.5,  the  SWR  computed  using  the  FDTD  method  is  compared  to  the 
measurements  and  theory  reported  in  [6.14].  The  FDTD  voltages  and  currents  were 
monitored  at  z  =  0.0  cm.  Circles  labeled  FDTD  voltage  and  FDTD  current  indicate  that 
the  SWR  was  computed  using  either  voltage  or  current.  The  results  from  the  FDTD 
method  appear  to  capture  the  general  trend  of  the  measured  response;  however,  the 
agreement  is  not  that  good.  This  is  likely  due  to  the  fact  that  several  defining  features  of 
the  feed  and  the  dipole  were  approximated. 

When  applying  this  method  to  the  balun-fed  folded  dipole  with  an  orthogonal 
ground  plane,  a  key  point  is  that  the  equivalent  surface  totally  encloses  the  radiating 
element.  This  presents  a  problem  since  the  mesh  used  for  calculating  radiation  patterns 
has  a  discretization  of  Ax  =  0.04  cm,  Ay  =  0.03  cm,  and  Az  =  0.05  cm.  Since  the 
orthogonal  ground  plane  has  physical  dimensions  of  40.64  cm  by  40.64  cm,  the  mesh 
would  have  to  be  more  than  1,016  cells  in  the  x-direction  and  more  than  1,355  cells  in  the 
y-direction.  These  numbers  are  each  an  order  of  magnitude  too  large  for  computers  with 
128  MB  of  RAM.  However,  by  truncating  the  mesh  in  the  x-  and  y-directions,  so  that  the 
orthogonal  ground  plane  extends  to  the  edge  of  the  mesh  in  the  x-  and  y-directions,  the 
ground  plane  appears  to  be  of  infinite  extent  as  far  as  the  mesh  is  concerned. 
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Figure  6.5  Standing-wave  ratio  of  the  balun-fed  folded  dipole.  FDTD  calculation  with 
reference  plane  at  z  =  -0.5  cm. 


Although  this  model  does  prevent  electric  fields  from  propagating  around  a  finite 
perfectly  conducting  plane,  it  will  not  lead  to  correct  radiation  patterns  in  the  lower  half 
plane  because  the  equivalent  surface  that  attempts  to  surround  the  radiator  does  not 
totally  enclose  it.  Instead,  the  orthogonal  ground  plane  intersects  the  equivalent  surface 
and  introduces  undesirable  effects  for  angles  near  grazing.  However,  it  does  provide  a 
useful  approximation  for  the  far  field  in  the  upper  half  plane. 

Radiation  patterns  comparing  the  co-polarized  and  cross-polarized  fields  with 
measurements  [6.14]  for  single  frequency  excitations  are  shown  in  Figs.  6.6  and  6.7.  In 
Fig.  6.6,  the  frequency  of  excitation  was  3.7  GHz;  there  is  good  agreement  between  the 
FDTD  results  and  the  measurements,  although  the  FDTD  cross-polarized  field  does  miss 
a  null.  Figure  6.7  compares  the  radiation  patterns  of  the  measurements  at  3.7  GHz  and 
the  FDTD  simulation  at  3.4  GHz.  The  agreement  between  the  FDTD  results  and 
measurements  is  better  at  this  frequency  than  in  Fig.  6.6.  Overall,  the  agreement  between 
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the  measurements  and  the  FDTD  results  is  very  good  considering  that  the  numerical 
model  has  a  much  smaller  orthogonal  ground  plane  than  the  structure  that  was  measured. 


(b) 

Figure  6.6  Radiation  patterns  at  frequency  3.7  GHz.  Co-polarized  field  is  labeled  e-co. 
Cross-polarized  field  is  labeled  e-cr.  (a)  E-plane,  (b)  H-plane. 
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(b) 

Figure  6.7  Radiation  patterns.  Measurements  are  at  frequency  3.7  GHz,  and  FDTD 
results  are  at  3.4  GHz.  Co-polarized  field  is  labeled  e-co.  Cross-polarized 
field  is  labeled  e-cr.  (a)  E-plane,  (b)  H-plane. 


6.4  Microstrip-Line  Fed  Patch  Antenna 

A  microstrip-line  fed  patch  antenna  was  built  and  measured.  S',,  and  a  radiation 
pattern  are  calculated  using  the  FDTD  model,  and  the  results  are  compared  with 
measurements.  The  antenna,  shown  in  Fig.  6.8,  is  center  fed  and  has  a  finite  ground 
plane.  The  presence  of  the  finite  ground  plane  makes  this  problem  well-suited  for  FDTD, 
whereas  moment  method  techniques  have  difficulty  handling  the  finite  ground  plane 
when  calculating  the  radiation  pattern.  Measurements  of  Sn  were  performed  on  the 
HP8510  network  analyzer.  TRL  standards,  shown  in  Fig.  6.9,  were  built  and  used  to 
calibrate  the  network  analyzer. 


Microstrip  Patch  Antenna 


Figure  6.8  Microstrip-line  fed  patch  antenna  with  finite  ground  plane. 
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LT  =  60.0  mm 

LR  =  30.0  mm 

Figure  6.9  TRL  standards  used  to  calibrate  network  analyzer. 


The  Prony  technique  was  used  to  calculate  5,,  at  the  reference  plane.  The 
magnitude  and  phase  of  Su  are  plotted  in  Figs.  6.10  and  6.11.  The  agreement  between 
the  measurements  and  the  FDTD  results  is  excellent: 

The  far  field  was  calculated  as  detailed  in  the  previous  section.  The  microstrip 
line  was  excited  at  the  resonance  frequency.  The  radiation  pattern  in  the  yz-plane  is 
shown  in  Fig.  6.12.  Measured  values  and  FDTD  results  are  both  normalized  to  30  dB  at 
0=0°  for  plotting  purposes.  Again,  the  agreement  between  the  results  is  excellent.  Note 
that  there  is  significant  power  in  the  z  <  0  half-space. 


Frequency  (GHz) 


Figure  6.10  Magnitude  of  5„  for  the  microstrip  patch  antenna. 


Frequency  (GHz) 


Figure  6.11  Phase  of  Sn  in  degrees  for  the  microstrip  patch  antenna. 
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Figure  6.12  Far-field  radiation  pattern  in  the  yz-plane.  lE-phil  vs.  theta. 

6.5  Coaxial-Line  Fed  Microstrip  Patch  Antenna 

Not  only  does  the  nonuniform  orthogonal  approach  result  in  considerable  savings 
of  memory  and  accelerate  computation  times,  it  also  allows  for  the  accurate  modeling  of 
rectangular  geometries.  In  this  section,  a  microstrip-line  fed  patch  antenna  is  analyzed  in 
order  to  demonstrate  the  reduction  in  required  memory.  Next,  coax  fed  patches  are 
analyzed,  demonstrating  the  flexibility  of  the  nonuniform  mesh. 

The  microstrip-line  fed  patch  antenna  investigated  in  this  section  is  similar  to  that 
shown  in  Fig.  6.8.  The  patch  is  4  cm  by  4  cm  and  the  feed  line  is  10  cm  long  and  0.4  cm 
wide.  The  substrate  is  0.15875  cm  thick  with  a  relative  dielectric  constant  of  2.55  and  the 
ground  plane  is  infinite.  This  geometry  was  modeled  using  a  uniform  mesh  that  was  120 
cells  by  180  cells  by  30  cells  with  Ax  =  Ay  =  0.1  cm  and  Az  =  0.0529166  cm,  and  the 
simulation  ran  for  8192  time  steps.  Using  the  nonuniform  feature  in  the  z-direction,  it 
was  possible  to  model  the  same  space  with  15  cells  as  opposed  to  30  cells.  Since  the 
number  of  unknowns  was  cut  in  half,  the  simulation  was  run  for  16384  time  steps  in  the 
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same  amount  of  computation  time  as  the  uniform  case.  Doubling  the  length  of  the  time 
signal  gives  twice  as  much  resolution  in  the  frequency  domain.  Figure  6.13  shows  the 
magnitude  of  Sn  computed  using  the  uniform  and  nonuniform  methods.  Note  that  the 
nonuniform  method  shows  a  deeper  null  that  is  missed  by  the  uniform  method  due  to  the 
shorter  time  signature.  Increased  resolution  in  the  frequency  domain  could  be  obtained 
with  the  same  total  computer  time  by  allowing  the  cells  to  grow  in  the  x-  and  y-directions 
as  well. 


Figure  6.13  Magnitude  of  Sn  of  microstrip-line  fed  patch  antenna.  Uniform  simulation 
had  8192  time  steps  and  nonuniform  simulation  had  16384  time  steps. 

It  is  straightforward  to  model  a  microstrip  line  fed  patch  antenna.  To  obtain  a 
50- 12  line,  a  circuit  simulator  is  used  to  determine  the  width  of  the  conductor  given  the 
substrate  thickness  and  dielectric  constant.  The  coax  fed  patch  is  slightly  more  involved. 
The  circular  inner  and  outer  conductors  are  approximated  by  square  conductors.  This 
approximation  is  reasonable  if  the  proper  impedance  can  be  obtained.  Choices  for  feed 
impedances  are  severely  limited  using  the  uniform  method.  Equation  (6.5)  gives  the  ratio 
of  the  widths  of  the  outer  and  inner  conductors  for  a  dielectric  filled  square  coaxial  line 
with  impedance  Z0  [6.15]. 
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w_  0.9259 

d  ~  gZ^/59.37 


(6.5) 


For  the  case  of  a  50-  Q  line  with  er  =  3.6,  d/w  =  5.3383.  Clearly,  this  would  be  difficult 
to  realize  using  a  small  number  of  cells  and  the  uniform  approach.  However,  with  the 
nonuniform  method,  this  ratio  can  be  obtained  using  five  cells  in  each  direction  as  shown 
in  Fig.  6. 14. 
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Case  2:  d/w  =  5.620 
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Case  4:  d/w  =  6.625 


Figure  6.14  Cross-sections  of  four  different  meshes  to  describe  a  50-  Q  square  coaxial 
line.  In  all  four  cases,  er  =  3.6. 
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Unfortunately,  the  FDTD  modeling  of  the  ratios  determined  by  (6.5)  does  not  lead 
to  the  desired  impedance,  so  the  ratio  of  the  widths  of  the  inner  and  outer  conductors  was 
varied  as  shown  in  Fig.  6.14,  and  the  resulting  impedances  are  shown  in  Fig.  6.15.  These 
results  were  obtained  using  the  Prony  technique  described  previously  to  extract  the 
reflection  due  to  the  imperfect  absorbing  boundaries. 


Figure  6.15  Magnitude  of  impedance  of  square  coaxial  line. 


A  coax  fed  patch  antenna  is  shown  in  Fig.  6.16.  The  dimensions  of  the  patch  and 
the  position  of  the  feed  point  would  require  a  prohibitively  large  computational  domain 
using  a  uniform  mesh.  With  the  nonuniform  mesh,  the  geometry  is  accurately  modeled 
with  a  reasonable  number  of  unknowns.  The  mesh  is  1 1 1  cells  by  1 1 1  cells  by  93  cells, 
and  there  is  approximately  60  mm  between  the  patch  antenna  and  the  absorbing 
boundaries. 
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Figure  6.16  Coax  fed  microstrip  patch  antenna.  The  dielectric  constant  is  3.6.  (a)  Top 
view,  (b)  Side  view. 


The  microstrip  patch  antenna  was  designed  for  circular  polarization,  built,  and 
measured.  Tuning  stubs  were  used  to  enhance  the  performance  of  the  antenna,  but 
records  of  the  sizes  and  locations  of  the  stubs  were  not  kept.  Measured  results  are  shown 
in  Fig.  6.17,  where  the  frequency  starts  at  f  j  =  1 .48  GHz  and  stops  at  f2  =  1.60  GHz,  and 
the  location  of  the  reference  plane  is  unknown.  This  antenna  was  modeled  using 
nonuniform  orthogonal  FDTD,  and  results  are  shown  in  Figs.  6.18  -  6.20,  where  the 


frequency  starts  at  fj  =  1.4420  GHz  and  stops  at  {2  =  1.5461  GHz,  and  the  reference  plane 
is  located  at  the  junction  between  the  coaxial  feed  and  the  ground  plane.  The  input 
impedance  for  the  antenna  is  shown  in  Fig.  6.18.  There  is  a  loop  near  resonance,  so  the 
difference  between  the  lengths  of  the  sides  has  to  be  reduced  for  better  performance.  A 
4.2  mm  by  2.7  mm  tuning  stub  was  placed  on  the  50.5  mm  side,  effectively  increasing  the 
length  of  the  shorter  side.  The  input  impedance  for  the  antenna  with  this  tuning  stub  is 
shown  in  Fig.  6.19.  The  cusp  in  the  input  impedance  is  not  very  deep,  suggesting  that  the 
tuning  stub  is  too  large.  The  stub  was  reduced  to  2.8  mm  by  2.7  mm.  The  input 
impedance  for  the  antenna  with  the  smaller  tuning  stub  is  shown  in  Fig.  6.20.  The 
difference  between  the  measured  and  calculated  resonant  frequencies  is  less  than  1.5%. 
Further  refinements  of  the  position  and  size  of  the  tuning  stub  would  further  improve  the 
agreement  between  the  measured  results  and  the  FDTD  results. 


Figure  6.17  Measured  input  impedance  of  coax  fed  patch  antenna.  Frequencies  fj  and  f2 
are  1.48  GHz  and  1.60  GHz,  respectively. 
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Figure  6.18  Input  impedance  of  coax  fed  patch  antenna.  Frequencies  fj  and  f2  are 
1.442  GHz  and  1.5461  GHz,  respectively. 


Figure  6.19  Input  impedance  of  coax  fed  patch  antenna  with  tuning  stub.  Frequencies  fj 
and  f2  are  1.442  GHz  and  1.5461  GHz,  respectively. 
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Si i  Tuning  Stub  2 


Figure  6.20  Input  impedance  of  coax  fed  patch  antenna  with  tuning  stub.  Frequencies  fj 
and  f2  are  1.442  GHz  and  1.5461  GHz,  respectively. 


6.6  Conclusions 

The  FDTD  method  was  applied  to  analyze  a  very  complicated  geometry,  namely  a 
printed  circuit  folded  dipole  with  an  integrated  balun.  The  numerical  results  from  the 
FDTD  model  were  compared  with  theoretical  and  measured  results.  The  results  for  the 
radiation  pattern  were  in  excellent  agreement  with  the  measured  data.  However,  the 
SWR  calculations  were  much  more  sensitive  to  the  modeling  approximations.  Therefore, 
a  center  fed  microstrip  line  fed  patch  antenna  was  constructed,  measured,  and  modeled 
with  FDTD.  The  FDTD  results  were  shown  to  be  in  excellent  agreement  with  the 
experimental  measurements.  The  nonuniform  orthogonal  method  was  used  in  order  to 
accurately  model  coaxially  fed  microstrip  patch  antennas.  The  nonuniformity  allows 
much  greater  accuracy  in  describing  the  geometry  as  compared  with  the  uniform  method, 
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which  requires  that  patch  dimensions  and  features  be  integral  multiples  of  the  cell 
discretization.  Such  a  restriction  is  entirely  impractical.  Although  the  nonuniform  grid 
gives  enhanced  modeling  capabilites,  there  is  room  for  improvement  in  the  modeling  and 
calculation  of  input  impedance  of  coaxially  fed  patch  antennas  using  the  FDTD  method, 
specifically  for  design  purposes. 
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CHAPTER  7 

CONCLUSIONS  AND  FUTURE  WORK 


The  nonuniform  orthogonal  FDTD  method  was  developed  and  demonstrated  on  a 
wide  variety  of  problems.  The  nonuniform  orthogonal  FDTD  method  is  superior  to  the 
conventional  FDTD  method  in  several  aspects.  The  nonuniform  discretization  provides 
great  flexibility  in  modeling  geometries,  whereas  dimensions  must  be  integral  multiples 
of  a  fixed  cell  discretization  in  the  conventional  method.  Moreover,  the  ability  of  the 
mesh  to  be  dense  in  areas  of  interest  and  then  gradually  expand  to  a  relatively  coarse 
mesh  makes  it  possible  to  efficiently  analyze  two  classes  of  problems.  These  classes  of 
problems  are  (a)  geometries  with  fine  features,  and  (b)  large  geometries.  Futhermore,  the 
orthogonality  of  the  method  preserves  the  speed  of  the  conventional  method. 

Chapter  2  presented  the  nonuniform  orthogonal  FDTD  method.  Update  equations 
were  derived  from  the  general  curvilinear  FDTD  update  equations  and  shown  to  be  of  the 
same  form  as  the  conventional  update  equations,  thus  demonstrating  the  preservation  of 
the  speed  of  the  update  equations.  Memory  requirements  were  shown  to  be  only 
marginally  larger  than  for  the  conventional  method.  This  is  significant  when  compared 
with  the  curvilinear  method,  which  is  expensive  in  terms  of  both  memory  requirements 
and  loss  of  speed.  The  error  associated  with  the  nonuniform  grid  was  discussed.  A 
numerical  example  demonstrated  that  gradual  growth  rates  in  the  mesh  discretization  lead 
to  errors  that  can  be  maintained  at  acceptably  low  levels.  A  more  rigorous  approach  to 
the  analysis  of  the  nonuniform  grid  error  is  suggested  for  future  work. 

A  dipole  radiating  in  the  presence  of  lossy  layered  media  was  analyzed,  and 
shown  to  give  excellent  results  when  compared  with  analytic  and  moment  method  results. 
The  nonuniform  orthogonal  FDTD  method  resulted  in  significant  savings  of  memory  and 
computation  time  when  compared  with  the  conventional  FDTD  method.  The  nonuniform 
orthogonal  FDTD  method  was  also  applied  to  the  finline  waveguide  structure.  The 


FDTD  results  obtained  were  shown  to  be  in  good  agreement  with  results  from  the  regular 
solution  of  the  singular  integral  equation. 
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Chapter  3  provided  a  numerical  analysis  of  the  dispersive  boundary  condition 
(DBC)  applied  to  nonuniform  grids.  Stability  and  dc  offsets  were  discussed,  and  the 
DBC  was.  tested  with  radiating  and  guided  wave  problems.  Error  analysis  was  carried  out 
in  the  time  domain,  the  frequency  domain,  and  spatially.  Analysis  of  Liao's  ABC  applied 
to  nonuniform  grids  is  a  suggested  direction  for  further  study. 

Chapter  4  presented  a  variety  of  methods  for  calculating  the  frequency-dependent 
characteristics  of  microstrip  lines,  striplines,  and  discontinuities.  A  Prony  technique  with 
the  FDTD  method  was  presented  and  applied  to  two-port  scattering  parameter 
calculations.  A  local  mesh  refinement  technique  was  introduced  for  triangular 
metallization  and  demonstrated  to  increase  the  accuracy  of  the  FDTD  method.  The 
FDTD  method  and  Prony  technique  were  combined  to  analyze  the  complicated  transition 
from  a  grounded  coplanar  waveguide  to  a  microstrip  line.  Further  analysis  of  the 
microstrip  line  to  the  grounded  coplanar  waveguide  transition  is  suggested.  Microstrip  or 
coplanar  waveguide  stubs  should  be  able  to  provide  a  good  match  at  frequencies  of 
interest. 

The  methods  of  Chapter  4  were  extended  in  Chapter  5  to  handle  a  broader  class  of 
problems.  The  new  method  was  presented,  followed  by  numerical  examples  analyzing 
coupled  symmetric  and  asymmetric  lines.  A  multidielectric  asymmetric  coupled  stripline 
problem  was  also  analyzed.  A  possible  avenue  for  research  involves  applying  the  general 
pencil-of-function  method  in  a  manner  similar  to  that  presented  in  Chapter  5.  This 
method  should  permit  the  extraction  of  general  modal  information  for  arbitrary  coupled 
line  structures. 

Chapter  6  discussed  the  application  of  the  FDTD  method  to  analyze  microstrip 
antennas.  A  near-field  to  far-field  transformation  was  discussed,  followed  by  the  analysis 
of  a  complicated  balun-fed  folded  dipole  antenna.  The  benefits  of  using  the  nonuniform 


137 


orthogonal  FDTD  method  were  demonstrated  through  the  analysis  of  microstrip  patch 
antennas.  The  results  obtained  were  compared  with  experimental  measurements. 
Excellent  results  were  obtained  for  the  microstrip  line  fed  patch  antenna.  However, 
improvements  can  be  made  to  calculate  the  input  impedance  more  accurately  when 
modeling  coaxially  fed  microstrip  patch  antennas. 
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